#include "context.h"
#include "evaluator_utils.h"
using namespace NTL;

constexpr uint64_t FOO_MAX_LEN = UINT64_MAX; // 定义最大长度常量
uint64_t Context::BASECONVERT = 0;
/*
 * @brief Construct a new Context:: Context object
 *
 * @param logN:多项式次数
 * @param logp: 单个模数q位数
 * @param L: 模数Q长度
 * @param K: 模数P长度
 * @param T: 模数T长度
 * @param gamma: 基转换分组数
 * @param h: 秘钥稀疏度
 * @param sigma: 高斯分布标准差
 */
Context::Context(long logN, long logp, long L, long K, long T, int gamma, int lambda, int dNum, long h, double sigma)
        : logN(logN), logp(logp), L(L), K(K), T(T), gamma(gamma), lambda(lambda), dNum(dNum), h(h), sigma(sigma)
{
    if(lambda == 0) {
        if ((logN == 15) && (L == 18)) {
            GetParameters(params::logN15QP891);
        }
        else if ((logN == 12) && (L == 28)) {
            GetParameters(params::logN12QP1725);
        }
        else if ((logN == 15) && (L == 28)) {
            GetParameters(params::logN15QP1725);
        }
        else if ((logN == 15) && (L == 15)) {
            GetParameters(params::logN15QP881);
        }
        else {
            cout << "没有满足要求的参数！！" << endl;
        }
        if (L == 0 || L > FOO_MAX_LEN) {
            throw std::domain_error("Invalid value for L.");
        }
        if (K == 0 || K > FOO_MAX_LEN) {
            throw std::domain_error("Invalid value for K.");
        }
        if (T == 0 || T > FOO_MAX_LEN) {
            throw std::domain_error("Invalid value for T.");
        }
    }
    else if ((lambda == 128) || (lambda == 192)){
        GetParametersLambda(logN, lambda, dNum);
    }
    else {
        throw std::domain_error("Invalid value for lambda.");
    }
    L = this->L;
    K = this->K;
    T = this->T;
    logp = this->logp;
    gamma = this->gamma;
    this->dNum = L / (K * 1.0);
    dNum = this->dNum;
    BASECONVERT = L; 
    N = 1L << logN;
    M = N << 1;
    logNh = logN - 1;
    Nh = N >> 1;
    this->p = 1L << logp;
    blockLen = ceil((L + K) / (gamma * 1.0));
    if (blockLen == 0 || blockLen > FOO_MAX_LEN) {
        throw std::runtime_error("Invalid value for blockLen.");
    }

    slots = N / 2;
    uVals = new complex<double>[slots]; // Context::Encode中的预分配内存，在Context::EncodePQ复用
    tmpCode = new uint64_t[N * L];          // Context::Decode中的预分配内存，Context::Back中复用
    tmp3 = new uint64_t[L << logN];     // Context::BackAndEqual，Context::Back，Context::RaiseAndEqual中复用
    tmpDecodeSingle = new uint64_t[N];  // Context::DecodeSingle的预分配内存
    raRaiseAndEqual = new uint64_t[(L + K) << logN]; // Context::RaiseAndEqual中的预分配内存
    raReScale = new uint64_t[(L + K) << logN];       // Context::ReScale中的预分配内存
    raBackAndEqual = new uint64_t[(L + K) << logN];  // Context::BackAndEqual中的预分配内存
    raReScaleAndEqual = new uint64_t[L << logN];     // Context::ReScaleAndEqual中的预分配内存
    tmpLeftRot = new uint64_t[L << logN];            // Context::LeftRot中的预分配内存

    KN = K * N;
    TN = T * N;
    ZZ tProduct = conv<ZZ>(1);
    for (int i = 0; i < T; i++) {
        tProduct *= conv<ZZ>(tVec[i]);
    }
    tProduct /= conv<ZZ>(2);

    halfTPolyPQ = new uint64_t[(L + K) << logN];
    polyQPool = new uint64_t[L << logN];
    polyTBootStrap = new uint64_t *[blockLen];
    for (int i = 0; i < blockLen; i++) {
        polyTBootStrap[i] = new uint64_t[T << logN];
    }

    // for Homomorphic Operator
    ex = new uint64_t[(L + K) << logN];
    sxSx = new uint64_t[(L + K) << logN];
    sxConj = new uint64_t[L << logN];
    new uint64_t[L << logN]();
    new uint64_t[L << logN]();

    new uint64_t[L << logN]();
    new uint64_t[L << logN]();
    new uint64_t[L << logN]();
    new uint64_t[L << logN]();
    new uint64_t[L << logN]();

    // Modulus
    qVecInv = new float[L]();
    twoqVec = new uint64_t[L]();
    qrVec = new __uint128_t[L];
    qrVecHigh = new uint64_t[L]();
    qrVecLow = new uint64_t[L]();
    qdVec = new uint64_t[L];
    qInvVec = new uint64_t[L];
    qRoots = new uint64_t[L];
    qRootsInv = new uint64_t[L];
    qRootPows = new uint64_t *[L];
    qRootPowsInv = new uint64_t *[L];
    NInvModq = new uint64_t[L];
    halfTPolyQi = new uint64_t[L];
    qRootPowsShoup = new uint64_t *[L]; 
    qRootPowsInvShoup = new uint64_t *[L]; 
    NInvModqShoup = new uint64_t[L]; 


    for (long i = 0; i < L; i++) {
        qVecInv[i] = 1 / qVec[i];

        halfTPolyQi[i] = conv<uint64_t>(tProduct % conv<ZZ>(qVec[i]));
        std::fill_n(halfTPolyPQ + ((K + i) << logN), N, halfTPolyQi[i]);
        // --------------------------Barrett-params------------------------------
        qrVec[i] = ~static_cast<__uint128_t>(0) / qVec[i];
        twoqVec[i] = static_cast<uint64_t>(2 * qVec[i]);
        qrVecHigh[i] = static_cast<uint64_t>(qrVec[i] >> 64);
        qrVecLow[i] = static_cast<uint64_t>(qrVec[i]);
        // --------------------------Barrett-params------------------------------
        qdVec[i] = qVec[i] << 1;
        qRoots[i] = FindMthRootOfUnity(M, qVec[i]);
        qRootsInv[i] = InvMod(qRoots[i], qVec[i]);
        NInvModq[i] = InvMod((uint64_t)N, qVec[i]);
        qInvVec[i] = Inv(qVec[i]);
        qRootPows[i] = new uint64_t[N];
        qRootPowsInv[i] = new uint64_t[N];
        uint64_t power = static_cast<uint64_t>(1);
        uint64_t powerInv = static_cast<uint64_t>(1);
        qRootPowsShoup[i] = new uint64_t[N]();          // shoup参数
        qRootPowsInvShoup[i] = new uint64_t[N]();        // shoup参数
        NInvModqShoup[i] = xShoup(NInvModq[i], qVec[i]); // shoup

        for (long j = 0; j < N; j++) {
            uint64_t jPrime = BitReverse(static_cast<uint32_t>(j)) >> (32 - logN);
            qRootPows[i][jPrime] = power;
            qRootPowsShoup[i][jPrime] = xShoup(power, qVec[i]); // shoup参数
            qRootPowsInv[i][jPrime] = powerInv;
            qRootPowsInvShoup[i][jPrime] = xShoup(powerInv, qVec[i]); // shoup参数

            if (j < N - 1) {
                MulMod(power, power, qRoots[i], qVec[i]);
                MulMod(powerInv, powerInv, qRootsInv[i], qVec[i]);
            }
        }
    }

    twopVec = new uint64_t[L]();
    prVecHigh = new uint64_t[L]();
    prVecLow = new uint64_t[L]();
    prVec = new __uint128_t[L];
    pdVec = new uint64_t[K];
    pInvVec = new uint64_t[K];
    pRoots = new uint64_t[K];
    pRootsInv = new uint64_t[K];
    pRootPows = new uint64_t *[K];
    pRootPowsInv = new uint64_t *[K];
    NInvModp = new uint64_t[K];
    NInvModpShoup = new uint64_t[K]; // shoup

    pRootPowsShoup = new uint64_t *[K];   // shoup参数
    pRootPowsInvShoup = new uint64_t *[K]; // shoup参数

    halfTPolyPi = new uint64_t[K];

    for (long i = 0; i < K; i++) {
        halfTPolyPi[i] = conv<uint64_t>(tProduct % conv<ZZ>(pVec[i]));
        std::fill_n(halfTPolyPQ + (i << logN), N, halfTPolyPi[i]);
        // --------------------------Barrett-params------------------------------
        prVec[i] = ~static_cast<__uint128_t>(0) / pVec[i];
        twopVec[i] = static_cast<uint64_t>(2 * pVec[i]);
        prVecHigh[i] = static_cast<uint64_t>(prVec[i] >> 64);
        prVecLow[i] = static_cast<uint64_t>(prVec[i]);
        // --------------------------Barrett-params------------------------------
        pdVec[i] = pVec[i] << 1;
        pRoots[i] = FindMthRootOfUnity(M, pVec[i]);
        pRootsInv[i] = InvMod(pRoots[i], pVec[i]);
        NInvModp[i] = InvMod((uint64_t)N, pVec[i]);
        pRootPows[i] = new uint64_t[N];
        pRootPowsInv[i] = new uint64_t[N];
        pInvVec[i] = Inv(pVec[i]);
        uint64_t power = static_cast<uint64_t>(1);
        uint64_t powerInv = static_cast<uint64_t>(1);
        pRootPowsShoup[i] = new uint64_t[N]();           // shoup参数
        pRootPowsInvShoup[i] = new uint64_t[N]();         // shoup参数
        NInvModpShoup[i] = xShoup(NInvModp[i], pVec[i]); // shoup
        for (long j = 0; j < N; j++) {
            uint64_t jPrime = BitReverse(static_cast<uint32_t>(j)) >> (32 - logN);
            pRootPows[i][jPrime] = power;
            pRootPowsShoup[i][jPrime] = xShoup(power, pVec[i]); // shoup参数
            pRootPowsInv[i][jPrime] = powerInv;
            pRootPowsInvShoup[i][jPrime] = xShoup(powerInv, pVec[i]); // shoup参数

            if (j < N - 1) {
                MulMod(power, power, pRoots[i], pVec[i]);
                MulMod(powerInv, powerInv, pRootsInv[i], pVec[i]);
            }
        }
    }

    qHatModq = new uint64_t *[L]; // [l][i] (phat_i)_l mod p_i
    for (long l = 0; l < L; l++) {
        qHatModq[l] = new uint64_t[l + 1];
        for (long i = 0; i < l + 1; i++) {
            qHatModq[l][i] = 1;
            for (long j = 0; j < i; j++) {
                uint64_t temp = qVec[j] % qVec[i];
                MulMod(qHatModq[l][i], qHatModq[l][i], temp, qVec[i]);
            }
            for (long j = i + 1; j < l + 1; j++) {
                uint64_t temp = qVec[j] % qVec[i];
                MulMod(qHatModq[l][i], qHatModq[l][i], temp, qVec[i]);
            }
        }
    }

    qHatInvModq = new uint64_t *[L]; // [l][i] (phat_i)_l^-1 mod p_i
    for (long l = 0; l < L; l++) {
        qHatInvModq[l] = new uint64_t[l + 1];
        for (long i = 0; i < l + 1; i++) {
            qHatInvModq[l][i] = InvMod(qHatModq[l][i], qVec[i]);
        }
    }

    pHatModp = new uint64_t[K]; // [k] qhat_k mod q_k

    for (long k = 0; k < K; k++) {
        pHatModp[k] = 1;
        for (long j = 0; j < k; j++) {
            uint64_t temp = pVec[j] % pVec[k];
            MulMod(pHatModp[k], pHatModp[k], temp, pVec[k]);
        }
        for (long j = k + 1; j < K; j++) {
            uint64_t temp = pVec[j] % pVec[k];
            MulMod(pHatModp[k], pHatModp[k], temp, pVec[k]);
        }
    }

    pHatInvModp = new uint64_t[K]; // [k] qhat_k^-1 mod q_k
    pHatInvModpShoup = new uint64_t[K];
    for (long k = 0; k < K; k++) {
        pHatInvModp[k] = InvMod(pHatModp[k], pVec[k]);
        pHatInvModpShoup[k] = xShoup(pHatInvModp[k], pVec[k]);
    }

    qHatModp = new uint64_t **[L]; // [l] [i] [k]  (phat_i)_l mod q_k
    for (long l = 0; l < L; l++) {
        qHatModp[l] = new uint64_t *[l + 1];
        for (long i = 0; i < l + 1; i++) {
            qHatModp[l][i] = new uint64_t[K];
            for (long k = 0; k < K; k++) {
                qHatModp[l][i][k] = 1;
                for (long j = 0; j < i; j++) {
                    uint64_t temp = qVec[j] % pVec[k];
                    MulMod(qHatModp[l][i][k], qHatModp[l][i][k], temp, pVec[k]);
                }
                for (long j = i + 1; j < l + 1; j++) {
                    uint64_t temp = qVec[j] % pVec[k];
                    MulMod(qHatModp[l][i][k], qHatModp[l][i][k], temp, pVec[k]);
                }

            }
        }
    }

    pHatModq = new uint64_t *[K]; // [k][i] qhat_k mod p_i
    for (long k = 0; k < K; k++) {
        pHatModq[k] = new uint64_t[L];
        for (long i = 0; i < L; i++) {
            pHatModq[k][i] = 1;
            for (long s = 0; s < k; s++) {
                uint64_t temp = pVec[s] % qVec[i];
                MulMod(pHatModq[k][i], pHatModq[k][i], temp, qVec[i]);
            }
            for (long s = k + 1; s < K; s++) {
                uint64_t temp = pVec[s] % qVec[i];
                MulMod(pHatModq[k][i], pHatModq[k][i], temp, qVec[i]);
            }
        }
    }

    PModq = new uint64_t[L]; // [i] qprod mod p_i
    for (long i = 0; i < L; i++) {
        PModq[i] = 1;
        for (long k = 0; k < K; k++) {
            uint64_t temp = pVec[k] % qVec[i];
            MulMod(PModq[i], PModq[i], temp, qVec[i]);
        }
    }

    PInvModq = new uint64_t[L]; // [i] qprod^-1 mod p_i
    PInvModqShoup = new uint64_t[L];
    for (long i = 0; i < L; i++) {
        PInvModq[i] = InvMod(PModq[i], qVec[i]);
        PInvModqShoup[i] = xShoup(PInvModq[i], qVec[i]);
    }

    QModp = new uint64_t *[L];
    for (long i = 0; i < L; i++) {
        QModp[i] = new uint64_t[K];
        for (long k = 0; k < K; k++) {
            QModp[i][k] = 1;
            for (long j = 0; j < i + 1; j++) {
                uint64_t temp = qVec[j] % pVec[k];
                MulMod(QModp[i][k], QModp[i][k], temp, pVec[k]);
            }
        }
    }
    QInvModp = new uint64_t *[L];
    for (long i = 0; i < L; i++) {
        QInvModp[i] = new uint64_t[K];
        for (long k = 0; k < K; k++) {
            QInvModp[i][k] = InvMod(QModp[i][k], pVec[k]);
        }
    }

    qInvModq = new uint64_t *[L]; // [i][j] p_i^-1 mod p_j
    for (long i = 0; i < L; i++) {
        qInvModq[i] = new uint64_t[L];
        for (long j = 0; j < i; j++) {
            qInvModq[i][j] = InvMod(qVec[i], qVec[j]);
        }
        for (long j = i + 1; j < L; j++) {
            qInvModq[i][j] = InvMod(qVec[i], qVec[j]);
        }
    }

    rotGroup = new long[Nh];
    long fivePows = 1;
    for (long i = 0; i < Nh; i++) {
        rotGroup[i] = fivePows;
        fivePows *= 5;
        fivePows %= M;
    }

    ksiPows = new complex<double>[M + 1];
    for (long j = 0; j < M; j++) {
        double angle = 2.0 * M_PI * j / M;
        ksiPows[j].real(cos(angle));
        ksiPows[j].imag(sin(angle));
    }

    ksiPows[M] = ksiPows[0];
    tVecInv = new float[T]();
    trVec = new __uint128_t[T];
    twotVec = new uint64_t[L]();
    trVecHigh = new uint64_t[L]();
    trVecLow = new uint64_t[L]();

    tInvVec = new uint64_t[T];

    tdVec = new uint64_t[T];
    tRoots = new uint64_t[T];

    tRootsInv = new uint64_t[T];
    tRootPows = new uint64_t *[T];
    tRootPowsInv = new uint64_t *[T];

    NInvModt = new uint64_t[T];


    halfTPolyti = new uint64_t[T];
    halfTPolyT = new uint64_t[T << logN];

    tRootPowsShoup = new uint64_t *[T];    // shoup
    tRootPowsInvShoup = new uint64_t *[T]; // shoup
    NInvModtShoup = new uint64_t[T];       // shoup


    for (long i = 0; i < T; i++) {
        tVecInv[i] = 1 / tVec[i];
        halfTPolyti[i] = conv<uint64_t>(tProduct % conv<ZZ>(tVec[i]));
        std::fill_n(halfTPolyT + (i << logN), N, halfTPolyti[i]);

        trVec[i] = ~static_cast<__uint128_t>(0) / tVec[i];

        twotVec[i] = static_cast<uint64_t>(2 * tVec[i]);
        trVecHigh[i] = static_cast<uint64_t>(trVec[i] >> 64);
        trVecLow[i] = static_cast<uint64_t>(trVec[i]);
        tdVec[i] = tVec[i] << 1;
        tRoots[i] = FindMthRootOfUnity(M, tVec[i]);

        tRootsInv[i] = InvMod(tRoots[i], tVec[i]);
        NInvModt[i] = InvMod((uint64_t)N, tVec[i]);

        tRootPows[i] = new uint64_t[N];

        tRootPowsInv[i] = new uint64_t[N];
        tInvVec[i] = Inv(tVec[i]);
        uint64_t power = static_cast<uint64_t>(1);
        uint64_t powerInv = static_cast<uint64_t>(1);
        tRootPowsShoup[i] = new uint64_t[N];             // shoup
        tRootPowsInvShoup[i] = new uint64_t[N];          // shoup
        NInvModtShoup[i] = xShoup(NInvModt[i], tVec[i]); // shoup

        for (long j = 0; j < N; j++) {
            uint64_t jprime = BitReverse(static_cast<uint32_t>(j)) >> (32 - logN);
            tRootPows[i][jprime] = power;
            tRootPowsShoup[i][jprime] = xShoup(power, tVec[i]); // shoup
            tRootPowsInv[i][jprime] = powerInv;
            tRootPowsInvShoup[i][jprime] = xShoup(powerInv, tVec[i]);

            if (j < N - 1) {
                MulMod(power, power, tRoots[i], tVec[i]);
                MulMod(powerInv, powerInv, tRootsInv[i], tVec[i]);
            }
        }
    }

    for (auto i = 0; i < K; i++) {
        PQ.emplace_back(pVec[i]);
        pqrVecHigh.emplace_back(prVecHigh[i]);
        pqrVecLow.emplace_back(prVecLow[i]);
    }
    for (auto i = 0; i < L; i++) {
        PQ.emplace_back(qVec[i]);
        pqrVecHigh.emplace_back(qrVecHigh[i]);
        pqrVecLow.emplace_back(qrVecLow[i]);

        QP.emplace_back(qVec[i]);
        twoqpVec.emplace_back(twoqVec[i]);
        qprVecHigh.emplace_back(qrVecHigh[i]);
        qprVecLow.emplace_back(qrVecLow[i]);
    }

    for (auto i = 0; i < K; i++) {
        QP.emplace_back(pVec[i]);
        twoqpVec.emplace_back(twopVec[i]);
        qprVecHigh.emplace_back(prVecHigh[i]);
        qprVecLow.emplace_back(prVecLow[i]);
    }

    if (gamma == 0) {
        throw std::runtime_error("Division by zero is not allowed.");
    }
    auto num_slices = (PQ.size() - 1) / gamma + 1;
    for (auto i = 0; i < blockLen; i++) {
        auto offset = i * gamma;
        auto length = gamma;
        if (i == blockLen - 1) {
            length = (L + K) - offset;
        }
        vector<uint64_t> tmp;
        for (auto j = 0; j < length; j++) {
            tmp.emplace_back(PQ[offset + j]);
        }
        PQSlices.emplace_back(tmp);
    }

    PQSlicesHatModPQ.resize(blockLen);
    for (auto i = 0; i < blockLen; i++) {
        vector<uint64_t> tmp;
        auto length = PQSlices[i].size();
        auto offset = i * gamma;
        PQSlicesHatModPQ[i].resize(length);
        for (auto j = 0; j < length; j++) {
            PQSlicesHatModPQ[i][j] = 1;
            uint64_t modVec = 0;
            if (offset + j < K) {
                modVec = pVec[offset + j];
            } else {
                modVec = qVec[offset + j - K];
            }
            for (auto m = 0; m < j; m++) {
                uint64_t divVec = 0;
                if (offset + m < K) {
                    divVec = pVec[offset + m];
                } else {
                    divVec = qVec[offset + m - K];
                }
                uint64_t temp = divVec % modVec;
                MulMod(PQSlicesHatModPQ[i][j], PQSlicesHatModPQ[i][j], temp, modVec);
            }
            for (auto m = j + 1; m < length; m++) {
                uint64_t divVec = 0;
                if (offset + m < K) {
                    divVec = pVec[offset + m];
                } else {
                    divVec = qVec[offset + m - K];
                }
                uint64_t temp = divVec % modVec;
                MulMod(PQSlicesHatModPQ[i][j], PQSlicesHatModPQ[i][j], temp, modVec);
            }
        }
    }

    PQSlicesHatInvModPQ.resize(blockLen);
    PQSlicesHatInvModPQShoup.resize(blockLen);
    for (auto i = 0; i < blockLen; i++) {
        auto length = PQSlices[i].size();
        auto offset = i * gamma;
        PQSlicesHatInvModPQ[i].resize(length);
        PQSlicesHatInvModPQShoup[i].resize(length);
        for (auto j = 0; j < length; j++) {
            uint64_t modVec = 0;
            if (offset + j < K) {
                modVec = pVec[offset + j];
            } else {
                modVec = qVec[offset + j - K];
            }
            PQSlicesHatInvModPQ[i][j] = InvMod(PQSlicesHatModPQ[i][j], modVec);
            PQSlicesHatInvModPQShoup[i][j] = xShoup(PQSlicesHatInvModPQ[i][j], modVec);
        }
    }

    PQSlicesHatModt.resize(blockLen);
    for (auto i = 0; i < blockLen; i++) {
        auto length = PQSlices[i].size();
        PQSlicesHatModt[i].resize(length);
        auto offset = i * gamma;
        for (auto j = 0; j < length; j++) {
            PQSlicesHatModt[i][j].resize(T);
            for (auto m = 0; m < T; m++) {
                PQSlicesHatModt[i][j][m] = 1;
                for (auto o = 0; o < j; o++) {
                    uint64_t divMod = 0;
                    if (offset + o < K) {
                        divMod = pVec[offset + o];
                    } else {
                        divMod = qVec[offset + o - K];
                    }
                    uint64_t temp = divMod % tVec[m];
                    MulMod(PQSlicesHatModt[i][j][m], PQSlicesHatModt[i][j][m], temp, tVec[m]);
                }
                for (auto o = j + 1; o < length; o++) {
                    uint64_t divMod = 0;
                    if (offset + o < K) {
                        divMod = pVec[offset + o];
                    } else {
                        divMod = qVec[offset + o - K];
                    }
                    uint64_t temp = divMod % tVec[m];
                    MulMod(PQSlicesHatModt[i][j][m], PQSlicesHatModt[i][j][m], temp, tVec[m]);
                }
            }
        }
    }

    PQSlicesModt.resize(blockLen);
    for (auto i = 0; i < blockLen; i++) {
        auto length = PQSlices[i].size();
        PQSlicesModt[i].resize(T);
        auto offset = i * gamma;
        for (auto j = 0; j < T; j++) {
            PQSlicesModt[i][j] = 1;
            for (auto m = 0; m < length; m++) {
                uint64_t curVec = 0;
                if (offset + m < K) {
                    curVec = pVec[offset + m];
                } else {
                    curVec = qVec[offset + m - K];
                }
                uint64_t temp = curVec % tVec[j];
                MulMod(PQSlicesModt[i][j], PQSlicesModt[i][j], temp, tVec[j]);
            }
        }
    }
    auto size = (L - 1) / K + 1;
    num_slices = (L - 1) / size + 1;
    qSlicesHatModq.resize(dNum);
    for (auto i = 0; i < dNum; i++) {
        auto offset = i * K;
        auto length = K;
        if (i == dNum - 1) {
            length = L - offset;
        }
        qSlicesHatModq[i].resize(length);
        for (auto j = 0; j < length; j++) {
            qSlicesHatModq[i][j].resize(j + 1);
            for (auto k = 0; k < j + 1; k++) {
                qSlicesHatModq[i][j][k] = 1;
                for (auto m = 0; m < k; m++) {
                    uint64_t temp = qVec[offset + m] % qVec[offset + k];
                    MulMod(qSlicesHatModq[i][j][k], qSlicesHatModq[i][j][k], temp, qVec[offset + k]);
                }
                for (auto m = k + 1; m < j + 1; m++) {
                    uint64_t temp = qVec[offset + m] % qVec[offset + k];
                    MulMod(qSlicesHatModq[i][j][k], qSlicesHatModq[i][j][k], temp, qVec[offset + k]);
                }
            }
        }
    }

    qSlicesHatInvModq.resize(dNum);
    qSlicesHatInvModqShoup.resize(dNum);
    for (auto i = 0; i < dNum; i++) {
        auto offset = i * K;
        auto length = K;
        if (i == dNum - 1) {
            length = L - offset;
        }
        qSlicesHatInvModq[i].resize(length);
        qSlicesHatInvModqShoup[i].resize(length);
        for (auto j = 0; j < length; j++) {
            qSlicesHatInvModq[i][j].resize(j + 1);
            qSlicesHatInvModqShoup[i][j].resize(j + 1);
            for (auto k = 0; k < j + 1; k++) {
                qSlicesHatInvModq[i][j][k] = InvMod(qSlicesHatModq[i][j][k], qVec[offset + k]);
                qSlicesHatInvModqShoup[i][j][k] = xShoup(qSlicesHatInvModq[i][j][k], qVec[offset + k]);
            }
        }
    }

    qSlicesHatModt.resize(dNum);
    for (auto i = 0; i < dNum; i++) {
        auto offset = i * K;
        auto length = K;
        if (i == dNum - 1) {
            length = L - offset;
        }
        qSlicesHatModt[i].resize(length);
        for (auto j = 0; j < length; j++) {
            qSlicesHatModt[i][j].resize(j + 1);
            for (auto k = 0; k < j + 1; k++) {
                qSlicesHatModt[i][j][k].resize(T);
                for (auto m = 0; m < T; m++) {
                    qSlicesHatModt[i][j][k][m] = 1;
                    for (auto n = 0; n < k; n++) {
                        uint64_t temp = qVec[offset + n] % tVec[m];
                        MulMod(qSlicesHatModt[i][j][k][m], qSlicesHatModt[i][j][k][m], temp, tVec[m]);
                    }
                    for (auto n = k + 1; n < j + 1; n++) {
                        uint64_t temp = qVec[offset + n] % tVec[m];
                        MulMod(qSlicesHatModt[i][j][k][m], qSlicesHatModt[i][j][k][m], temp, tVec[m]);
                    }
                }
            }
        }
    }

    QSlicesModt.resize(dNum);
    for (auto i = 0; i < dNum; i++) {
        auto offset = i * K;
        auto length = K;
        if (i == dNum - 1) {
            length = L - offset;
        }
        QSlicesModt[i].resize(length);
        for (auto j = 0; j < length; j++) {
            QSlicesModt[i][j].resize(T);
            for (auto k = 0; k < T; k++) {
                QSlicesModt[i][j][k] = 1;
                for (auto m = 0; m < j + 1; m++) {
                    uint64_t temp = qVec[offset + m] % tVec[k];
                    MulMod(QSlicesModt[i][j][k], QSlicesModt[i][j][k], temp, tVec[k]);
                }
            }
        }
    }

    tHatModt.resize(T);
    for (auto k = 0; k < T; k++) {
        tHatModt[k] = 1;
        for (auto j = 0; j < k; j++) {
            uint64_t temp = tVec[j] % tVec[k];
            MulMod(tHatModt[k], tHatModt[k], temp, tVec[k]);
        }
        for (auto j = k + 1; j < T; j++) {
            uint64_t temp = tVec[j] % tVec[k];
            MulMod(tHatModt[k], tHatModt[k], temp, tVec[k]);
        }
    }

    tHatInvModt.resize(T);
    tHatInvModtShoup.resize(T);
    for (auto k = 0; k < T; k++) {
        tHatInvModt[k] = InvMod(tHatModt[k], tVec[k]);
        tHatInvModtShoup[k] = xShoup(tHatInvModt[k], tVec[k]);
    }

    tHatModPQ.resize(T);
    for (auto k = 0; k < T; k++) {
        tHatModPQ[k].resize(L + K);
        for (auto i = 0; i < K; i++) {
            tHatModPQ[k][i] = 1;
            for (auto s = 0; s < k; s++) {
                uint64_t temp = tVec[s] % pVec[i];
                MulMod(tHatModPQ[k][i], tHatModPQ[k][i], temp, pVec[i]);
            }
            for (auto s = k + 1; s < T; s++) {
                uint64_t temp = tVec[s] % pVec[i];
                MulMod(tHatModPQ[k][i], tHatModPQ[k][i], temp, pVec[i]);
            }
        }
        for (auto i = K; i < K + L; i++) {
            tHatModPQ[k][i] = 1;
            for (auto s = 0; s < k; s++) {
                uint64_t temp = tVec[s] % qVec[i - K];
                MulMod(tHatModPQ[k][i], tHatModPQ[k][i], temp, qVec[i - K]);
            }
            for (auto s = k + 1; s < T; s++) {
                uint64_t temp = tVec[s] % qVec[i - K];
                MulMod(tHatModPQ[k][i], tHatModPQ[k][i], temp, qVec[i - K]);
            }
        }
    }

    TModPQ.resize(K + L);
    for (auto i = 0; i < K; i++) {
        TModPQ[i] = 1;
        for (auto k = 0; k < T; k++) {
            uint64_t temp = tVec[k] % pVec[i];
            MulMod(TModPQ[i], TModPQ[i], temp, pVec[i]);
        }
    }
    for (auto i = K; i < L + K; i++) {
        TModPQ[i] = 1;
        for (auto k = 0; k < T; k++) {
            uint64_t temp = tVec[k] % qVec[i - K];
            MulMod(TModPQ[i], TModPQ[i], temp, qVec[i - K]);
        }
    }

    // usually the bound of m_std * M is used, where M = 12 .. 40
    // we use M = 12 here, which corresponds to the probability of roughly 2^(-100)
    constexpr double acc{ 5e-32 };
    //    constexpr double M{ noexcept(sqrt(-2 * log(acc))) };
    int fin{ static_cast<int>(ceil(sigma * M)) };
    
    mVals.clear();
    mVals.reserve(fin);
    double variance{ 2 * sigma * sigma };
    double cusum{ 0.0 };
    for (int x = 1; x <= fin; ++x) {
        cusum += exp(-(static_cast<double>(x * x) / variance));
        mVals.push_back(cusum);
    }
    double m_a = 1.0 / (2 * cusum + 1.0);
    #pragma omp parallel for
    for (int x = 0; x < fin; ++x)
        mVals[x] *= m_a;
}

/*
 * @brief ：Construct a new Context:: Context object
 * 
 * @param logN ：多项式次数
 * @param lambda ：安全等级
 * @param d ：模数Q与模数P的比值
 * @return Context ：预计算参数
 */
Context Context::InitContext(long logN, int lambda, int d)
{
    long L = 15;
    long K = 1;
    long T = 5;
    long logp = 55;
    long gamma = 4;

    return Context(logN, logp, L, K, T, gamma, lambda, d);
}

Context Context::InitContextBootStrapping(long logN)
{
    long L = 28;
    long K = 4;
    long T = 8;
    long logp = 53;
    long gamma = 4;
    return Context(logN, logp, L, K, T, gamma);
}


/*
 * @brief 对复数数组进行位反转操作
 *
 * @param vals: 输入的复数数组，数组中的每个元素在反转过程中会被重新排列。
 * @param size: 数组的大小
 */
void Context::ArrayBitReverse(complex<double> *vals, const long size)
{
    for (long i = 1, j = 0; i < size; i++) {
        long bit = size >> 1;
        for (; j >= bit; bit >>= 1) {
            j -= bit;
        }
        j += bit;
        if (i < j) {
            swap(vals[i], vals[j]);
        }
    }
}

/*
 * @brief 对无符号 64 位整数数组进行位反转操作
 *
 * @param vals: 输入的无符号 64 位整数数组，数组中的每个元素在反转过程中会被重新排列。
 * @param size: 数组的大小
 */
void Context::ArrayBitReverse(uint64_t *vals, const long size)
{
    for (long i = 1, j = 0; i < size; i++) {
        long bit = size >> 1;
        for (; j >= bit; bit >>= 1) {
            j -= bit;
        }
        j += bit;
        if (i < j) {
            swap(vals[i], vals[j]);
        }
    }
}


/*
 * @brief 快速傅里叶变换（FFT）算法

 * @param vals: 输入的复数数组，经过 FFT 变换后将被更新为傅里叶变换的结果
 * @param size: 数组的大小
 */
void Context::FFTSpecial(complex<double> *vals, const long size)
{
    ArrayBitReverse(vals, size);
    for (long len = 2; len <= size; len <<= 1) {
        for (long i = 0; i < size; i += len) {
            long lenH = len >> 1;
            long lenQ = len << 2;
            for (long j = 0; j < lenH; j++) {
                long idx = ((rotGroup[j] % lenQ)) * M / lenQ;
                complex<double> u = vals[i + j];
                complex<double> v = vals[i + j + lenH];
                v *= ksiPows[idx];
                vals[i + j] = u + v;
                vals[i + j + lenH] = u - v;
            }
        }
    }
}

/*
 * @brief：逆快速傅里叶变换（IFFT）算法
 *
 * @param vals：输入的复数数组
 * @param size：数组的大小
 */
void Context::FFTSpecialInvLazy(complex<double> *vals, const long size)
{
    for (long len = size; len >= 1; len >>= 1) {
        for (long i = 0; i < size; i += len) {
            long lenH = len >> 1;
            long lenQ = len << 2;
            for (long j = 0; j < lenH; j++) {
                long idx = (lenQ - (rotGroup[j] % lenQ)) * M / lenQ;
                complex<double> u = vals[i + j] + vals[i + j + lenH];
                complex<double> v = vals[i + j] - vals[i + j + lenH];
                v *= ksiPows[idx];
                vals[i + j] = u;
                vals[i + j + lenH] = v;
            }
        }
    }
    ArrayBitReverse(vals, size);
}

/*
 * @brief ：执行完整的逆快速傅里叶变换（IFFT），包含归一化处理。
 *
 * @param vals ：输入的复数数组
 * @param size ：数组的大小
 */
void Context::FFTSpecialInv(complex<double> *vals, const long size)
{
    if (size == 0) {
        throw std::runtime_error("Division by zero is not allowed.");
    }
    FFTSpecialInvLazy(vals, size);
    for (long i = 0; i < size; i++) {
        vals[i] /= size;
    }
}

/*
 * @brief ：对复数数组进行编码，将数据映射到多项式系数表示中。
 *
 * @param a ：输出的多项式数组，每个槽的复数值将被编码为整数并存储在该数组中。
 * @param v ：输入的复数数组，包含待编码的数据，长度为 slotsEncode。
 * @param slotsEncode ：表示数据的槽数，决定输入数组的大小。
 * @param l ：编码过程中使用的模数个数，用于多项式系数的映射。
 */
void Context::Encode(uint64_t *a, complex<double> *v, long slotsEncode, long l)
{
    if (slotsEncode == 0) {
        throw std::runtime_error("Division by zero is not allowed.");
    }
    std::copy(v, v + slotsEncode, uVals);

    long gap = Nh / slotsEncode;

    FFTSpecialInv(uVals, slotsEncode);

    for (long j = 0; j < slotsEncode; j++) {
        uVals[j] *= p;
    }
    {
        for (long i = 0; i < l; i++) {
            uint64_t *mi = a + i * N;
            for (long j = 0, jdx = Nh, idx = 0; j < slotsEncode; j++, jdx += gap, idx += gap) {
                double mir_value = uVals[j].real();
                double mii_value = uVals[j].imag();
                NTL::ZZ mir_zz, mii_zz;
                uint64_t mirr, miii;
                if (mir_value >= 0) {
                    mirr = static_cast<uint64_t>(mir_value);
                    mirr = mirr % qVec[i];
                    mi[idx] = mirr;
                } else {
                    mirr = static_cast<uint64_t>(-mir_value);
                    mirr = mirr % qVec[i];
                    mi[idx] = qVec[i] - mirr;
                }
                if (mii_value >= 0) {
                    miii = static_cast<uint64_t>(mii_value);
                    miii = miii % qVec[i];
                    mi[jdx] = miii;
                } else {
                    miii = static_cast<uint64_t>(-mii_value);
                    miii = miii % qVec[i];
                    mi[jdx] = qVec[i] - miii;
                }
            }
            QiNTTAndEqual(mi, i);
        }
    }
}

/*
 * @brief ：对复数数组进行多模数编码，将数据映射到环PQ多项式系数表示中。
 *
 * @param a： 输出的多项式数组，编码后的整数数据将存储到该数组中。
 * @param v： 输入的复数数组，包含待编码的数据，长度为 slotsEncode。
 * @param scale： 缩放因子，用于对输入数据进行放大。
 * @param slotsEncode： 数据槽的数量，决定输入数组的大小。
 * @param l： 编码中涉及的模数Q数量，用于多项式系数的映射。
 * @param k： 编码中涉及的模数P数量，用于多项式系数的映射。
 */
void Context::EncodePQ(uint64_t *a, complex<double> *v, double scale, long slotsEncode, long l, long k)
{
    if (slotsEncode == 0) {
        throw std::runtime_error("Division by zero is not allowed.");
    }
    std::copy(v, v + slotsEncode, uVals);

    long gap = Nh / slotsEncode;

    FFTSpecialInv(uVals, slotsEncode);

    for (long j = 0; j < slotsEncode; j++) {
        uVals[j] *= scale;
    }
    if (abs(uVals[0]) < p) {
        for (long i = 0; i < k; i++) {
            uint64_t *mi = a + i * N;
            for (long j = 0, jdx = Nh, idx = 0; j < slotsEncode; j++, jdx += gap, idx += gap) {
                long mir = uVals[j].real();
                long mii = uVals[j].imag();
                mi[idx] = mir >= 0 ? (uint64_t)mir : (uint64_t)(pVec[i] + mir);
                mi[jdx] = mii >= 0 ? (uint64_t)mii : (uint64_t)(pVec[i] + mii);
            }
            PiNTTAndEqual(mi, i);
        }
        for (long i = 0; i < l; i++) {
            uint64_t *mi = a + k * N + i * N;
            for (long j = 0, jdx = Nh, idx = 0; j < slotsEncode; j++, jdx += gap, idx += gap) {
                long mir = uVals[j].real();
                long mii = uVals[j].imag();
                mi[idx] = mir >= 0 ? (uint64_t)mir : (uint64_t)(qVec[i] + mir);
                mi[jdx] = mii >= 0 ? (uint64_t)mii : (uint64_t)(qVec[i] + mii);
            }
            QiNTTAndEqual(mi, i);
        }
    } else {
        for (long i = 0; i < k; i++) {
            uint64_t *mi = a + i * N;
            for (long j = 0, jdx = Nh, idx = 0; j < slotsEncode; j++, jdx += gap, idx += gap) {
                double mir_value = uVals[j].real();
                double mii_value = uVals[j].imag();
                NTL::ZZ mir_zz, mii_zz;
                uint64_t mirr, miii;
                if (mir_value >= 0) {
                    mirr = static_cast<uint64_t>(mir_value);
                    mirr = mirr % qVec[i];
                    mi[idx] = mirr;
                } else {
                    mirr = static_cast<uint64_t>(-mir_value);
                    mirr = mirr % qVec[i];
                    mi[idx] = qVec[i] - mirr;
                }
                if (mii_value >= 0) {
                    miii = static_cast<uint64_t>(mii_value);
                    miii = miii % qVec[i];
                    mi[jdx] = miii;
                } else {
                    miii = static_cast<uint64_t>(-mii_value);
                    miii = miii % qVec[i];
                    mi[jdx] = qVec[i] - miii;
                }
            }
            PiNTTAndEqual(mi, i);
        }
        for (long i = 0; i < l; i++) {
            uint64_t *mi = a + k * N + i * N;
            for (long j = 0, jdx = Nh, idx = 0; j < slotsEncode; j++, jdx += gap, idx += gap) {
                double mir_value = uVals[j].real();
                double mii_value = uVals[j].imag();
                NTL::ZZ mir_zz, mii_zz;
                uint64_t mirr, miii;
                if (mir_value >= 0) {
                    mirr = static_cast<uint64_t>(mir_value);
                    mirr = mirr % pVec[i];
                    mi[idx] = mirr;
                } else {
                    mirr = static_cast<uint64_t>(-mir_value);
                    mirr = mirr % pVec[i];
                    mi[idx] = pVec[i] - mirr;
                }
                if (mii_value >= 0) {
                    miii = static_cast<uint64_t>(mii_value);
                    miii = miii % pVec[i];
                    mi[jdx] = miii;
                } else {
                    miii = static_cast<uint64_t>(-mii_value);
                    miii = miii % pVec[i];
                    mi[jdx] = pVec[i] - miii;
                }
            }
            QiNTTAndEqual(mi, i);
        }
    }
}

void Context::Encode(uint64_t *ax, double *vals, long slots, long l) {}

/*
 * @brief ：对单个复数值进行编码，将其映射到整数多项式的系数表示中
 *
 * @param axEncode  ：输出多项式
 * @param val  ：输入的复数值
 * @param l  ：编码过程中使用的模数个数
 */
void Context::EncodeSingle(uint64_t *axEncode, complex<double> &val, long l)
{
    long long vr = val.real() * p;
    long long vi = val.imag() * p;

    for (long i = 0; i < l; i++) {
        uint64_t *ai = axEncode + i * N;
        while (vr < 0) {
            vr += qVec[i];
        }
        ai[0] = vr;

        while (vi < 0) {
            vi += qVec[i];
        }
        ai[Nh] = vi;
        QiNTTAndEqual(ai, i);
    }
}

/*
 * @brief ：对单个实数值进行编码
 *
 * @param ax ：输出多项式
 * @param val ：输入的实数值
 * @param l ：编码过程中使用的模数个数
 */
void Context::EncodeSingle(uint64_t *ax, double val, long l)
{
    complex<double> complexVal(val, 0.0);
    EncodeSingle(ax, complexVal, l);
}

/*
 * @brief ：对整数多项式进行解码，将其转换回复数形式的数据表示。
 *
 * @param a： 输入的多项式数组，包含经过编码和变换的数据。
 * @param v： 输出的复数数组，解码后的复数值将存储在此数组中。
 * @param slotsDecode： 数据槽的数量，表示复数数组的有效长度。
 * @param l： 多项式编码过程中使用的模数个数。
 * @param scale： 解码时的缩放因子，用于恢复原始复数值。
 */
void Context::Decode(uint64_t *a, complex<double> *v, long slotsDecode, long l, __uint128_t scale)
{
    if ((slotsDecode == 0) || (scale == 0)) {
        throw std::runtime_error("Division by zero is not allowed.");
    }
    std::copy(a, a + N * l, tmpCode);
    long gap = Nh / slotsDecode;

    INTTAndEqual(tmpCode, l);

    std::vector<ZZ> coeffsBigint(N);
    fill(coeffsBigint.begin(), coeffsBigint.end(), 0);

    PolyToBigintLvl(l, tmpCode, int(gap), coeffsBigint);

    ZZ Q_modules(1);
    for (int i = 0; i < l; i++) {
        Q_modules *= qVec[i];
    }
    ZZ Q_half = Q_modules / 2;
    int sign;
    ZZ c;
    ZZ ci;
    for (long j = 0, jdx = Nh, idx = 0; j < slotsDecode; j++, jdx += gap, idx += gap) {
        c = coeffsBigint[idx];
        ci = coeffsBigint[jdx];
        c = coeffsBigint[idx] < Q_half ? coeffsBigint[idx] : (coeffsBigint[idx] - Q_modules);
        ci = coeffsBigint[jdx] < Q_half ? coeffsBigint[jdx] : (coeffsBigint[jdx] - Q_modules);
        v[j].real(conv<double>(c) / scale);
        v[j].imag(conv<double>(ci) / scale);
    }
    FFTSpecial(v, slotsDecode);

}

/*
 * @brief ：将多项式的模数表示转换为大整数表示
 *
 * @param level ：模数的个数
 * @param p1 ：输入的多项式数组
 * @param gap ：每个数据槽之间的间隔
 * @param coeffsBigint ：输出的大整数向量
 */
void Context::PolyToBigintLvl(int level, uint64_t *p1, int gap, std::vector<ZZ> &coeffsBigint)
{
    uint64_t qi;

    std::vector<ZZ> crtReconstruction(level);

    ZZ QiB;
    ZZ tmp;
    ZZ modulusBigint(1);
    ZZ ModulusBigint = conv<ZZ>(1);
    for (int i = 0; i < level; i++) {
        ModulusBigint *= conv<ZZ>(qVec[i]);
    }
    for (int i = 0; i < level; i++) {
        qi = qVec[i];
        QiB = conv<ZZ>(qi);
        modulusBigint *= QiB;
        crtReconstruction[i] = ModulusBigint / QiB;
        tmp = crtReconstruction[i] % QiB;
        InvMod(tmp, tmp, QiB);
        tmp %= QiB;
        crtReconstruction[i] *= tmp;
    }

    for (int i = 0, j = 0; j < N; i++, j += gap) {
        tmp = 0;
        for (int k = 0; k < level; k++) {
            coeffsBigint[i] += conv<ZZ>(p1[k * N + j]) * crtReconstruction[k];
        }
        coeffsBigint[i] %= modulusBigint;
    }
}

/*
 * @brief ：对单个整数多项式进行解码，将其转换为复数形式的数值表示。
 *
 * @param ax ：输入的多项式数组
 * @param val ：输出的复数值
 * @param l ：模数个数
 */
void Context::DecodeSingle(uint64_t *ax, complex<double> &val)
{
    std::copy(ax, ax + N, tmpDecodeSingle);
    QiINTTAndEqual(tmpDecodeSingle, 0);

    uint64_t pr = qVec[0];
    uint64_t pr_2 = qVec[0] / 2;

    double vr = tmpDecodeSingle[0] <= pr_2 ?
                static_cast<double>(tmpDecodeSingle[0]) / p :
                (static_cast<double>(tmpDecodeSingle[0]) - static_cast<double>(pr)) / static_cast<double>(p);
    double vi = tmpDecodeSingle[Nh] <= pr_2 ?
                static_cast<double>(tmpDecodeSingle[Nh]) / p :
                (static_cast<double>(tmpDecodeSingle[Nh]) - static_cast<double>(pr)) / static_cast<double>(p);

    val.real(vr);
    val.imag(vi);
}

/*
 * @brief ：对单个整数多项式进行解码，将其转换为实数形式的数值表示。
 *
 * @param ax ：输入的多项式数组
 * @param val ：输出的实数值
 * @param l ：模数个数
 */
void Context::DecodeSingle(uint64_t *ax, double val)
{
    std::copy(ax, ax + N, tmpCode);
    QiINTTAndEqual(tmpCode, 0);
    uint64_t pr = qVec[0];
    uint64_t pr_2 = qVec[0] / 2;

    double vr = tmpCode[0] <= pr_2 ? static_cast<double>(tmpCode[0]) / p :
                static_cast<double>(tmpCode[0]) - static_cast<double>(pr) / static_cast<double>(p);
    val = vr;
}

/*
 * @brief ：对多项式在模数PQ下进行数论变换 (NTT)，并将结果存储到输出数组中。
 *
 * @param res ：输出数组
 * @param a ：输入的多项式数组
 * @param l： 模数Q层级数量。
 * @param k： 模数P层级数量。
 */
void Context::NTT(uint64_t *res, uint64_t *a, long l, long k)
{
    for (long index = 0; index < l; index++) {
        uint64_t *ai = a + (index << logN);
        uint64_t *resi = res + (index << logN);
        QiNTT(resi, ai, index);
    }

    for (long index = l; index < l + k; index++) {
        uint64_t *ai = a + (index << logN);
        uint64_t *resi = res + (index << logN);
        PiNTT(resi, ai, index - l);
    }
}


/*
 * @brief ：对多项式在模数Q下进行数论变换 (NTT)，并将结果存储到输出数组中。
 *
 * @param res ：输出数组
 * @param a ：输入的多项式数组
 * @param index ：模数层级的索引
 */
void Context::QiNTT(uint64_t *res, uint64_t *a, long index)
{
    std::copy(a, a + N, res);
    QiNTTAndEqual(res, index);
}

/*
 * @brief ：对多项式在模数P下进行数论变换 (NTT)，并将结果存储到输出数组中。
 *
 * @param res ：输出数组
 * @param a ：输入的多项式数组
 * @param index ：模数层级的索引
 */
void Context::PiNTT(uint64_t *res, uint64_t *a, long index)
{
    std::copy(a, a + N, res);
    PiNTTAndEqual(res, index);
}

/*
 * @brief ：在模数Q指定模数下，对多项式数组执行数论变换 (NTT)，并将结果存储回原数组中
 *
 * @param a ：输入/输出数组，包含待变换的多项式数据，结果会直接更新到该数组中。
 * @param index ：模数索引
 */
void Context::QiNTTAndEqual(uint64_t *a, long index)
{
    long t = N;
    long logt1 = logN;
    uint64_t q = qVec[index];
    int gap = N / 8;
    for (long m = 1; m < gap; m <<= 1) {
        t >>= 1;
        for (long i = 0; i < m; i++) {
            long j1 = i << logt1;
            long j2 = j1 + t;
            uint64_t W = qRootPows[index][m + i];
            uint64_t wShoup = qRootPowsShoup[index][m + i];
            uint64x2_t q2 = vdupq_n_u64(q);
            for (long j = j1; j < j2; j += 4) {
                uint64_t V[4];
                for (int k = 0; k < 4; k++) {
                    uint128_t vMul = static_cast<uint128_t>(a[j + t + k]) * wShoup;
                    uint64_t hi = vMul >> 64;
                    V[k] = a[j + t + k] * W - hi * q;
                    if (V[k] >= q) {
                        V[k] -= q;
                    }
                }
                uint64x2_t vVec[2];
                vVec[0] = vcombine_u64(vdup_n_u64(V[0]), vdup_n_u64(V[1]));
                vVec[1] = vcombine_u64(vdup_n_u64(V[2]), vdup_n_u64(V[3]));
                for (int k = 0; k < 2; k++) {
                    uint64x2_t ajtVec = vld1q_u64(&a[j + t + 2 * k]);
                    uint64x2_t ajVec = vld1q_u64(&a[j + 2 * k]);
                    uint64x2_t ajPlusQMinusV = vsubq_u64(vaddq_u64(ajVec, q2), vVec[k]);
                    uint64x2_t lessMask = vcltq_u64(ajVec, vVec[k]);
                    uint64x2_t ajtNew = vbslq_u64(lessMask, ajPlusQMinusV, vsubq_u64(ajVec, vVec[k]));
                    vst1q_u64(&a[j + t + 2 * k], ajtNew);

                    ajVec = vaddq_u64(ajVec, vVec[k]);
                    uint64x2_t gtMask = vcgeq_u64(ajVec, q2);
                    ajVec = vbslq_u64(gtMask, vsubq_u64(ajVec, q2), ajVec);
                    vst1q_u64(&a[j + 2 * k], ajVec);
                }
            }
        }
        logt1 -= 1;
    }

    for (long m = gap; m < N; m <<= 1) {
        t >>= 1;
        for (long i = 0; i < m; i++) {
            long j1 = i << logt1;
            long j2 = j1 + t;
            uint64_t W = qRootPows[index][m + i];
            uint64_t wShoup = qRootPowsShoup[index][m + i];
            for (long j = j1; j < j2; j++) {
                uint128_t v_w_mul = static_cast<uint128_t>(a[j + t]) * wShoup;
                uint64_t hi = v_w_mul >> 64;
                uint64_t V = a[j + t] * W - hi * q;
                if (V >= q) {
                    V -= q;
                }
                a[j + t] = a[j] < V ? a[j] + q - V : a[j] - V;
                a[j] += V;
                if (a[j] >= q) {
                    a[j] -= q;
                }
            }
        }
        logt1 -= 1;
    }
}


/*
 * @brief ：在模数P指定模数下，对多项式数组执行数论变换 (NTT)，并将结果存储回原数组中
 *
 * @param a ：输入/输出数组，包含待变换的多项式数据，结果会直接更新到该数组中。
 * @param index ：模数索引
 */
void Context::PiNTTAndEqual(uint64_t *a, long index)
{
    long t = N;
    long logt1 = logN;
    uint64_t pModule = pVec[index];
    int gap = N / 8;
    for (long m = 1; m < gap; m <<= 1) {
        t >>= 1;
        for (long i = 0; i < m; i++) {
            long j1 = i << logt1;
            long j2 = j1 + t;
            uint64_t W = pRootPows[index][m + i];
            uint64_t wShoup = pRootPowsShoup[index][m + i];
            uint64x2_t q2 = vdupq_n_u64(pModule);
            for (long j = j1; j < j2; j += 4) {
                uint128_t vMul[4];
                uint64_t hi[4];
                uint64_t V[4];
                for (int k = 0; k < 4; k++) {
                    vMul[k] = static_cast<uint128_t>(a[j + t + k]) * wShoup;
                    hi[k] = vMul[k] >> 64;
                    V[k] = a[j + t + k] * W - hi[k] * pModule;
                    if (V[k] >= pModule) {
                        V[k] -= pModule;
                    }
                }
                uint64x2_t vVec1 = vcombine_u64(vdup_n_u64(V[0]), vdup_n_u64(V[1]));
                uint64x2_t vVec2 = vcombine_u64(vdup_n_u64(V[2]), vdup_n_u64(V[3]));
                for (int k = 0; k < 2; k++) {
                    uint64x2_t ajtVec = vld1q_u64(&a[j + t + 2 * k]);
                    uint64x2_t ajVec = vld1q_u64(&a[j + 2 * k]);
                    uint64x2_t ajPlusQMinusV = vsubq_u64(vaddq_u64(ajVec, q2), k == 0 ? vVec1 : vVec2);
                    uint64x2_t lessMask = vcltq_u64(ajVec, k == 0 ? vVec1 : vVec2);
                    uint64x2_t ajtNew = vbslq_u64(lessMask, ajPlusQMinusV, vsubq_u64(ajVec, k == 0 ? vVec1 : vVec2));
                    vst1q_u64(&a[j + t + 2 * k], ajtNew);

                    ajVec = vaddq_u64(ajVec, k == 0 ? vVec1 : vVec2);
                    uint64x2_t gtMask = vcgeq_u64(ajVec, q2);
                    ajVec = vbslq_u64(gtMask, vsubq_u64(ajVec, q2), ajVec);
                    vst1q_u64(&a[j + 2 * k], ajVec);
                }
            }
        }
        logt1 -= 1;
    }

    for (long m = gap; m < N; m <<= 1) {
        t >>= 1;
        for (long i = 0; i < m; i++) {
            long j1 = i << logt1;
            long j2 = j1 + t;
            uint64_t W = pRootPows[index][m + i];
            uint64_t wShoup = pRootPowsShoup[index][m + i];
            for (long j = j1; j < j2; j++) {
                uint128_t v_w_mul = static_cast<uint128_t>(a[j + t]) * wShoup;
                uint64_t hi = v_w_mul >> 64;
                uint64_t V = a[j + t] * W - hi * pModule;
                if (V >= pModule) {
                    V -= pModule;
                }
                a[j + t] = a[j] < V ? a[j] + pModule - V : a[j] - V;
                a[j] += V;
                if (a[j] >= pModule) {
                    a[j] -= pModule;
                }
            }
        }
        logt1 -= 1;
    }
}

/*
 * @brief ：在模数T指定模数下，对多项式数组执行数论变换 (NTT)，并将结果存储回原数组中
 *
 * @param a ：输入/输出数组，包含待变换的多项式数据，结果会直接更新到该数组中。
 * @param index ：模数索引
 */
void Context::TiNTTAndEqual(uint64_t *a, long index)
{
    long t = N;
    long logt1 = logN;
    uint64_t ti = tVec[index];
    int gap = N / 8;
    for (long m = 1; m < gap; m <<= 1) {
        t >>= 1;
        for (long i = 0; i < m; i++) {
            long j1 = i << logt1;
            long j2 = j1 + t;
            uint64_t W = tRootPows[index][m + i];
            uint64_t wShoup = tRootPowsShoup[index][m + i];
            uint64x2_t q2 = vdupq_n_u64(ti);
            for (long j = j1; j < j2; j += 4) {
                uint128_t vMul[4];
                uint64_t hi[4];
                uint64_t V[4];
                for (int k = 0; k < 4; k++) {
                    vMul[k] = static_cast<uint128_t>(a[j + t + k]) * wShoup;
                    hi[k] = vMul[k] >> 64;
                    V[k] = a[j + t + k] * W - hi[k] * ti;
                    if (V[k] >= ti) {
                        V[k] -= ti;
                    }
                }
                uint64x2_t vVec1 = vcombine_u64(vdup_n_u64(V[0]), vdup_n_u64(V[1]));
                uint64x2_t vVec2 = vcombine_u64(vdup_n_u64(V[2]), vdup_n_u64(V[3]));
                for (int k = 0; k < 2; k++) {
                    uint64x2_t ajtVec = vld1q_u64(&a[j + t + 2 * k]);
                    uint64x2_t ajVec = vld1q_u64(&a[j + 2 * k]);
                    uint64x2_t ajPlusQMinusV = vsubq_u64(vaddq_u64(ajVec, q2), k == 0 ? vVec1 : vVec2);
                    uint64x2_t lessMask = vcltq_u64(ajVec, k == 0 ? vVec1 : vVec2);
                    uint64x2_t ajtNew = vbslq_u64(lessMask, ajPlusQMinusV, vsubq_u64(ajVec, k == 0 ? vVec1 : vVec2));
                    vst1q_u64(&a[j + t + 2 * k], ajtNew);

                    ajVec = vaddq_u64(ajVec, k == 0 ? vVec1 : vVec2);
                    uint64x2_t gtMask = vcgeq_u64(ajVec, q2);
                    ajVec = vbslq_u64(gtMask, vsubq_u64(ajVec, q2), ajVec);
                    vst1q_u64(&a[j + 2 * k], ajVec);
                }
            }
        }
        logt1 -= 1;
    }

    for (long m = gap; m < N; m <<= 1) {
        t >>= 1;
        for (long i = 0; i < m; i++) {
            long j1 = i << logt1;
            long j2 = j1 + t;
            uint64_t W = tRootPows[index][m + i];
            uint64_t wShoup = tRootPowsShoup[index][m + i];
            for (long j = j1; j < j2; j++) {
                uint128_t v_w_mul = static_cast<uint128_t>(a[j + t]) * wShoup;
                uint64_t hi = v_w_mul >> 64;
                uint64_t V = a[j + t] * W - hi * ti;
                if (V >= ti) {
                    V -= ti;
                }
                a[j + t] = a[j] < V ? a[j] + ti - V : a[j] - V;
                a[j] += V;
                if (a[j] >= ti) {
                    a[j] -= ti;
                }
            }
        }
        logt1 -= 1;
    }
}


/*
 * @brief ：对输入多项式数组在模数PQ下执行原地数论变换 (NTT)，并更新数组数据。
 *
 * @param a ：输入/输出数组
 * @param l： 模数Q层级数量。
 * @param k： 模数P层级数量。
 */
void Context::NTTAndEqual(uint64_t *a, long l, long k)
{
#pragma omp parallel for
    for (long index = 0; index < l; index++) {
        uint64_t *ai = a + (index << logN);
        QiNTTAndEqual(ai, index);
    }

    for (long index = l; index < l + k; index++) {
        uint64_t *ai = a + (index << logN);
        PiNTTAndEqual(ai, index - l);
    }
}

/*
 * @brief ：对输入多项式数组载模数T下执行原地数论变换 (NTT)，并更新数组数据。
 *
 * @param a ：输入/输出数组
 * @param t ：模数T层级数量。
 */
void Context::TNTTAndEqual(uint64_t *a, long t)
{
// #pragma omp parallel for
    for (long index = 0; index < t; index++) {
        uint64_t *ai = a + (index << logN);
        TiNTTAndEqual(ai, index);
    }
}

/*
 * @brief ：对多项式在模数Q下进行逆数论变换 (INTT)，并将结果存储到输出数组中。
 *
 * @param res ：输出数组
 * @param a ：输入的多项式数组
 * @param index ：模数层级的索引
 */
void Context::QiINTT(uint64_t *res, uint64_t *a, long index)
{
    std::copy(a, a + N, res);
    QiINTTAndEqual(res, index);
}

/*
 * @brief ：对多项式在模数P下进行逆数论变换 (INTT)，并将结果存储到输出数组中。
 *
 * @param res ：输出数组
 * @param a ：输入的多项式数组
 * @param index ：模数层级的索引
 */
void Context::PiINTT(uint64_t *res, uint64_t *a, long index)
{
    std::copy(a, a + N, res);
    PiINTTAndEqual(res, index);
}

/*
 * @brief ：在模数T指定模数下，对多项式数组执行逆数论变换 (INTT)，并将结果存储回原数组中
 *
 * @param a ：输入/输出数组，包含待变换的多项式数据，结果会直接更新到该数组中。
 * @param index ：模数索引
 */
void Context::TiINTTAndEqual(uint64_t *a, long index)
{
    uint64_t ti = tVec[index];
    uint64_t td = tdVec[index];
    uint64_t tInv = tInvVec[index];
    long t = 1;
    for (long m = N; m > 1; m >>= 1) {
        long j1 = 0;
        long halfm = m >> 1;
        for (long i = 0; i < halfm; i++) {
            long j2 = j1 + t - 1;
            uint64_t W = tRootPowsInv[index][halfm + i];
            uint64_t WShoup = tRootPowsInvShoup[index][halfm + i];

            for (long j = j1; j <= j2; j++) {
                uint64_t Ti = a[j] + td;
                Ti -= a[j + t];
                a[j] += a[j + t];
                if (a[j] >= td)
                    a[j] -= td;
                uint128_t vWMul = static_cast<uint128_t>(Ti) * WShoup;
                uint64_t hi = vWMul >> 64;
                uint64_t V = Ti * W - hi * ti;
                if (V >= ti)
                    V -= ti;
                a[j + t] = V;
            }
            j1 += (t << 1);
        }
        t <<= 1;
    }
    uint64_t NScale = NInvModt[index];
    uint64_t NScaleShoup = NInvModtShoup[index];
    uint64x2_t tiVec = vdupq_n_u64(ti);
    uint64x2_t NScaleVec = vdupq_n_u64(NScale);
    for (long i = 0; i < N; i += 2) {
        uint64x2_t aVec = vld1q_u64(&a[i]);
        uint64x2_t TiVec = vsubq_u64(aVec, tiVec);
        uint64x2_t mask = vcgeq_u64(aVec, tiVec);
        TiVec = vbslq_u64(mask, TiVec, aVec);
        uint64_t Ti1 = vgetq_lane_u64(TiVec, 0);
        uint128_t vWMul1 = static_cast<uint128_t>(Ti1) * NScaleShoup;
        uint64_t hi1 = vWMul1 >> 64;
        uint64_t V1 = Ti1 * NScale - hi1 * ti;
        if (V1 >= ti) {
            V1 -= ti;
        }
        uint64_t Ti2 = vgetq_lane_u64(TiVec, 1);
        uint128_t vWMul2 = static_cast<uint128_t>(Ti2) * NScaleShoup;
        uint64_t hi2 = vWMul2 >> 64;
        uint64_t V2 = Ti2 * NScale - hi2 * ti;
        if (V2 >= ti) {
            V2 -= ti;
        }
        uint64x2_t VVec = { V1, V2 };
        vst1q_u64(&a[i], VVec);
    }
}

/*
 * @brief ：在模数Q指定模数下，对多项式数组执行逆数论变换 (INTT)，并将结果存储回原数组中
 *
 * @param a ：输入/输出数组，包含待变换的多项式数据，结果会直接更新到该数组中。
 * @param index ：模数索引
 */
void Context::QiINTTAndEqual(uint64_t *a, long index)
{
    uint64_t q = qVec[index];
    uint64_t qd = qdVec[index];
    uint64_t qInv = qInvVec[index];
    long t = 1;
    for (long m = N; m > 1; m >>= 1) {
        long j1 = 0;
        long halfm = m >> 1;
        for (long i = 0; i < halfm; i++) {
            long j2 = j1 + t - 1;
            uint64_t W = qRootPowsInv[index][halfm + i];
            uint64_t WShoup = qRootPowsInvShoup[index][halfm + i];

            for (long j = j1; j <= j2; j++) {
                uint64_t Ti = a[j] + qd;
                Ti -= a[j + t];
                a[j] += a[j + t];
                if (a[j] >= qd)
                    a[j] -= qd;
                uint128_t vWMul = static_cast<uint128_t>(Ti) * WShoup;
                uint64_t hi = vWMul >> 64;
                uint64_t V = Ti * W - hi * q;
                if (V >= q)
                    V -= q;
                a[j + t] = V;
            }
            j1 += (t << 1);
        }
        t <<= 1;
    }
    uint64_t NScale = NInvModq[index];
    uint64_t NScaleShoup = NInvModqShoup[index];
    uint64x2_t qiVec = vdupq_n_u64(q);
    uint64x2_t NScaleVec = vdupq_n_u64(NScale);
    for (long i = 0; i < N; i += 2) {
        uint64x2_t aVec = vld1q_u64(&a[i]);
        uint64x2_t TiVec = vsubq_u64(aVec, qiVec);
        uint64x2_t mask = vcgeq_u64(aVec, qiVec);
        TiVec = vbslq_u64(mask, TiVec, aVec);
        uint64_t Ti1 = vgetq_lane_u64(TiVec, 0);
        uint128_t vWMul1 = static_cast<uint128_t>(Ti1) * NScaleShoup;
        uint64_t hi1 = vWMul1 >> 64;
        uint64_t V1 = Ti1 * NScale - hi1 * q;
        if (V1 >= q) {
            V1 -= q;
        }
        uint64_t Ti2 = vgetq_lane_u64(TiVec, 1);
        uint128_t vWMul2 = static_cast<uint128_t>(Ti2) * NScaleShoup;
        uint64_t hi2 = vWMul2 >> 64;
        uint64_t V2 = Ti2 * NScale - hi2 * q;
        if (V2 >= q) {
            V2 -= q;
        }
        uint64x2_t VVec = { V1, V2 };
        vst1q_u64(&a[i], VVec);
    }
}

/*
 * @brief ：在模数P指定模数下，对多项式数组执行逆数论变换 (INTT)，并将结果存储回原数组中
 *
 * @param a ：输入/输出数组，包含待变换的多项式数据，结果会直接更新到该数组中。
 * @param index ：模数索引
 */
void Context::PiINTTAndEqual(uint64_t *a, long index)
{
    uint64_t pModule = pVec[index];
    uint64_t pd = pdVec[index];
    uint64_t pInv = pInvVec[index];
    long t = 1;
    for (long m = N; m > 1; m >>= 1) {
        long j1 = 0;
        long halfm = m >> 1;
        for (long i = 0; i < halfm; i++) {
            long j2 = j1 + t - 1;
            uint64_t W = pRootPowsInv[index][halfm + i];
            uint64_t WShoup = pRootPowsInvShoup[index][halfm + i];

            for (long j = j1; j <= j2; j++) {
                uint64_t Ti = a[j] + pd;
                Ti -= a[j + t];
                a[j] += a[j + t];
                if (a[j] >= pd)
                    a[j] -= pd;
                uint128_t vWMul = static_cast<uint128_t>(Ti) * WShoup;
                uint64_t hi = vWMul >> 64;
                uint64_t V = Ti * W - hi * pModule;
                if (V >= pModule)
                    V -= pModule;
                a[j + t] = V;
            }
            j1 += (t << 1);
        }
        t <<= 1;
    }
    uint64_t NScale = NInvModp[index];
    uint64_t NScaleShoup = NInvModpShoup[index];

    uint64x2_t piVec = vdupq_n_u64(pModule);

    for (long i = 0; i < N; i += 2) {
        uint64x2_t aVec = vld1q_u64(&a[i]);
        uint64x2_t TiVec = vsubq_u64(aVec, piVec);
        uint64x2_t mask = vcgeq_u64(aVec, piVec);
        TiVec = vbslq_u64(mask, TiVec, aVec);
        uint64_t Ti1 = vgetq_lane_u64(TiVec, 0);
        uint128_t vWMul1 = static_cast<uint128_t>(Ti1) * NScaleShoup;
        uint64_t hi1 = vWMul1 >> 64;
        uint64_t V1 = Ti1 * NScale - hi1 * pModule;
        if (V1 >= pModule) {
            V1 -= pModule;
        }
        uint64_t Ti2 = vgetq_lane_u64(TiVec, 1);
        uint128_t vWMul2 = static_cast<uint128_t>(Ti2) * NScaleShoup;
        uint64_t hi2 = vWMul2 >> 64;
        uint64_t V2 = Ti2 * NScale - hi2 * pModule;
        if (V2 >= pModule) {
            V2 -= pModule;
        }
        uint64x2_t VVec = { V1, V2 };
        vst1q_u64(&a[i], VVec);
    }
}

/*
 * @brief ： 对输入多项式数组执行原地逆数论变换 (INTT)，并更新数组数据。
 *
 * @param a ：输入/输出数组
 * @param l： 模数Q层级数量。
 * @param k： 模数P层级数量。
 */
void Context::INTTAndEqual(uint64_t *a, long l, long k)
{
#pragma omp parallel for
    for (long index = 0; index < l; index++) {
        uint64_t *ai = a + (index << logN);
        QiINTTAndEqual(ai, index);
    }

    for (long index = l; index < l + k; index++) {
        uint64_t *ai = a + (index << logN);
        PiINTTAndEqual(ai, index - l);
    }
}

/*
 * @brief ：对输入多项式数组载模数T下执行原地逆数论变换 (INTT)，并更新数组数据。
 *
 * @param a ：输入/输出数组
 * @param t ：模数T层级数量。
 */
void Context::TINTTAndEqual(uint64_t *a, long t)
{
#pragma omp parallel for
    for (long index = 0; index < t; index++) {
        uint64_t *ai = a + (index << logN);
        TiINTTAndEqual(ai, index);
    }
}

/*
 * @brief ：对多项式在模数PQ下进行逆数论变换 (NTT)，并将结果存储到输出数组中。
 *
 * @param res ：输出数组
 * @param a ：输入的多项式数组
 * @param l： 模数Q层级数量。
 * @param k： 模数P层级数量。
 */
void Context::INTT(uint64_t *res, uint64_t *a, long l, long k)
{
#pragma omp parallel for
    for (long index = 0; index < l; index++) {
        uint64_t *ai = a + (index << logN);
        uint64_t *resi = res + (index << logN);
        QiINTT(resi, ai, index);
    }

    for (long index = l; index < l + k; index++) {
        uint64_t *ai = a + (index << logN);
        uint64_t *resi = res + (index << logN);
        PiINTT(resi, ai, index - l);
    }
}

/*
 * @brief 在模数Q指定模数下，对输入数组中的多项式系数进行取负操作。
 *
 * @param res： 输出数组，用于存储取负后的结果。
 * @param a： 输入数组，包含待取负的多项式系数。
 * @param index： 模数索引，指定当前取负操作所使用的模数。
 */
void Context::QiNegate(uint64_t *res, uint64_t *a, long index)
{
    for (long i = 0; i < N; i++) {
        res[i] = (a[i] == 0) ? 0 : qVec[index] - a[i];
    }
}

/*
 * @brief 在模数P指定模数下，对输入数组中的多项式系数进行取负操作。
 *
 * @param res： 输出数组，用于存储取负后的结果。
 * @param a： 输入数组，包含待取负的多项式系数。
 * @param index： 模数索引，指定当前取负操作所使用的模数。
 */
void Context::PiNegate(uint64_t *res, uint64_t *a, long index)
{
    for (long i = 0; i < N; i++) {
        res[i] = (a[i] == 0) ? 0 : pVec[index] - a[i];
    }
}

/*
 * @brief 在指定模数下，对输入多项式数组执行取负操作，并将结果存储到输出数组中。
 *
 * @param res： 输出数组，用于存储取负后的多项式系数。
 * @param a： 输入数组，包含待取负的多项式系数。
 * @param l： 处理的模数层级数量，前 `l` 层调用 `QiNegate` 进行取负操作。
 * @param k： 额外层级数量，后 `k` 层调用 `PiNegate` 进行取负操作。
 */
void Context::Negate(uint64_t *res, uint64_t *a, long l, long k)
{
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *resi = res + (i << logN);
        QiNegate(resi, ai, i);
    }

    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *resi = res + (i << logN);
        PiNegate(resi, ai, i - l);
    }
}

/*
 * @brief 在模数Q指定模数下，对输入数组中的多项式系数进行取负操作。
 *
 * @param a： 输入/输出数组。
 * @param index： 模数索引，指定当前取负操作所使用的模数。
 */
void Context::QiNegateAndEqual(uint64_t *a, long index)
{
    for (long i = 0; i < N; i++) {
        a[i] = (a[i] == 0) ? 0 : qVec[index] - a[i];
    }
}

/*
 * @brief 在模数P指定模数下，对输入数组中的多项式系数进行取负操作。
 *
 * @param a： 输入/输出数组。
 * @param index： 模数索引，指定当前取负操作所使用的模数。
 */
void Context::PiNegateAndEqual(uint64_t *a, long index)
{
    for (long i = 0; i < N; i++) {
        a[i] = (a[i] == 0) ? 0 : pVec[index] - a[i];
    }
}

/*
 * @brief 在指定模数下，对输入多项式数组执行取负操作，并将结果存储到输出数组中。
 *
 * @param a： 输入/输出数组。
 * @param l： 处理的模数层级数量，前 `l` 层调用 `QiNegate` 进行取负操作。
 * @param k： 额外层级数量，后 `k` 层调用 `PiNegate` 进行取负操作。
 */
void Context::NegateAndEqual(uint64_t *a, long l, long k)
{
    for (long i = 0; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        #pragma omp parallel for
        for (long n = 0; n < N; n++) {
            ai[n] = (ai[n] == 0) ? 0 : QP[i] - a[n];
        }
    }
}

/*
 * @brief 在模数Q指定模数下，对输入多项式数组的每个元素加上一个常数，并存储结果。
 *
 * @param res： 输出数组，用于存储加上常数后的多项式系数。
 * @param a： 输入数组，包含待加常数的多项式系数。
 * @param c： 要加上的常数值。
 * @param index： 模数索引，用于从 `qVec` 中获取当前模数。
 */
void Context::QiAddConst(uint64_t *res, uint64_t *a, uint64_t c, long index)
{
    for (long i = 0; i < N; i++) {
        AddMod(res[i], a[i], c, qVec[index]);
    }
}

/*
 * @brief 在模数P指定模数下，对输入多项式数组的每个元素加上一个常数，并存储结果。
 *
 * @param res： 输出数组，用于存储加上常数后的多项式系数。
 * @param a： 输入数组，包含待加常数的多项式系数。
 * @param c： 要加上的常数值。
 * @param index： 模数索引，用于从 `qVec` 中获取当前模数。
 */
void Context::PiAddConst(uint64_t *res, uint64_t *a, uint64_t c, long index)
{
    for (long i = 0; i < N; i++) {
        AddMod(res[i], a[i], c, pVec[index]);
    }
}

/*
 * @brief 在指定模数域下，对输入多项式数组的每个元素加上一个常数，并存储结果。
 *
 * @param res： 输出数组，用于存储加常数后的多项式系数。
 * @param a： 输入数组，包含待加常数的多项式系数。
 * @param c： 要加上的常数值。
 * @param l： 处理的模数层级数量，前 `l` 层调用 `QiAddConst` 进行加法操作。
 * @param k： 额外层级数量，后 `k` 层调用 `PiAddConst` 进行加法操作。
 */
void Context::AddConst(uint64_t *res, uint64_t *a, uint64_t c, long l, long k)
{
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *resi = res + (i << logN);
        QiAddConst(resi, ai, c, i);
    }

    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *resi = res + (i << logN);
        PiAddConst(resi, ai, c, i - l);
    }
}

/*
 * @brief 在指定模数下，对输入多项式数组的每个元素加上一个常数，结果直接更新到原数组中。
 *
 * @param a： 输入/输出数组，包含待加常数的多项式系数，结果将直接更新到此数组。
 * @param c： 要加上的常数值。
 * @param index： 模数索引，用于从 `qVec` 中获取当前的模数。
 */
void Context::QiAddConstAndEqual(uint64_t *a, uint64_t c, long index)
{
    #pragma omp parallel for
    for (long i = 0; i < N; i++) {
        AddMod(a[i], a[i], c, qVec[index]);
    }
}

/*
 * @brief 在指定额外模数下，对输入多项式数组的每个元素加上一个常数，结果直接更新到原数组中。
 *
 * @param a： 输入/输出数组，包含待加常数的多项式系数，结果将直接更新到此数组。
 * @param c： 要加上的常数值。
 * @param index： 额外模数索引，用于从 `pVec` 中获取当前的模数。
 */
void Context::PiAddConstAndEqual(uint64_t *a, uint64_t c, long index)
{
    for (long i = 0; i < N; i++) {
        AddMod(a[i], a[i], c, pVec[index]);
    }
}

void Context::AddConstAndEqual(uint64_t *a, uint64_t c, long l, long k)
{
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        QiAddConstAndEqual(ai, c, i);
    }

    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        PiAddConstAndEqual(ai, c, i - l);
    }
}

void Context::QiSubConst(uint64_t *res, uint64_t *a, uint64_t c, long index)
{
    for (long i = 0; i < N; i++) {
        SubMod(res[i], a[i], c, qVec[index]);
    }
}

void Context::PiSubConst(uint64_t *res, uint64_t *a, uint64_t c, long index)
{
    for (long i = 0; i < N; i++) {
        SubMod(res[i], a[i], c, pVec[index]);
    }
}

void Context::SubConst(uint64_t *res, uint64_t *a, uint64_t c, long l, long k)
{
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *resi = res + (i << logN);
        QiSubConst(resi, ai, c, i);
    }

    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *resi = res + (i << logN);
        PiSubConst(resi, ai, c, i - l);
    }
}

void Context::QiSubConstAndEqual(uint64_t *a, uint64_t c, long index)
{
    for (long i = 0; i < N; i++) {
        SubMod(a[i], a[i], c, qVec[index]);
    }
}

void Context::PiSubConstAndEqual(uint64_t *a, uint64_t c, long index)
{
    for (long i = 0; i < N; i++) {
        SubMod(a[i], a[i], c, pVec[index]);
    }
}

void Context::SubConstAndEqual(uint64_t *a, uint64_t c, long l, long k)
{
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        QiSubConstAndEqual(ai, c, i);
    }

    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        PiSubConstAndEqual(ai, c, i - l);
    }
}

void Context::QiAdd(uint64_t *res, uint64_t *a, uint64_t *b, long index)
{
    uint64x2_t q_vec = vdupq_n_u64(qVec[index]);

    for (long i = 0; i < N; i += 2) {
        uint64x2_t a_vec = vld1q_u64(a + i);
        uint64x2_t b_vec = vld1q_u64(b + i);

        uint64x2_t sum = vaddq_u64(a_vec, b_vec);
        uint64x2_t cmp = vcgtq_u64(sum, q_vec);
        uint64x2_t res_vec = vsubq_u64(sum, vandq_u64(cmp, q_vec));
        vst1q_u64(res + i, res_vec);
    }
}


void Context::PiAdd(uint64_t *res, uint64_t *a, uint64_t *b, long index)
{
    uint64x2_t p_vec = vdupq_n_u64(pVec[index]);
#pragma omp parallel for
    for (long i = 0; i < N; i += 2) {
        uint64x2_t a_vec = vld1q_u64(a + i);
        uint64x2_t b_vec = vld1q_u64(b + i);

        uint64x2_t sum = vaddq_u64(a_vec, b_vec);
        uint64x2_t cmp = vcgtq_u64(sum, p_vec);
        uint64x2_t res_vec = vsubq_u64(sum, vandq_u64(cmp, p_vec));
        vst1q_u64(res + i, res_vec);
    }
}

void Context::Add(uint64_t *res, uint64_t *a, uint64_t *b, long l, long k)
{
#pragma omp parallel for
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        uint64_t *resi = res + (i << logN);
        QiAdd(resi, ai, bi, i);
    }

    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        uint64_t *resi = res + (i << logN);
        PiAdd(resi, ai, bi, i - l);
    }
}


void Context::QiAddAndEqual(uint64_t *a, uint64_t *b, long index)
{
    uint64x2_t q_vec = vdupq_n_u64(qVec[index]);

    for (long i = 0; i < N; i += 2) {
        uint64x2_t a_vec = vld1q_u64(a + i);
        uint64x2_t b_vec = vld1q_u64(b + i);
        uint64x2_t sum = vaddq_u64(a_vec, b_vec);
        uint64x2_t cmp = vcgtq_u64(sum, q_vec);
        uint64x2_t res_vec = vsubq_u64(sum, vandq_u64(cmp, q_vec));
        vst1q_u64(a + i, res_vec);
    }
}


void Context::PiAddAndEqual(uint64_t *a, uint64_t *b, long index)
{
    uint64x2_t p_vec = vdupq_n_u64(pVec[index]);

    for (long i = 0; i < N; i += 2) {
        uint64x2_t a_vec = vld1q_u64(a + i);
        uint64x2_t b_vec = vld1q_u64(b + i);
        uint64x2_t sum = vaddq_u64(a_vec, b_vec);
        uint64x2_t cmp = vcgtq_u64(sum, p_vec);
        uint64x2_t res_vec = vsubq_u64(sum, vandq_u64(cmp, p_vec));
        vst1q_u64(a + i, res_vec);
    }
}

void Context::TiAddAndEqual(uint64_t *a, uint64_t *b, long index)
{
    uint64x2_t t_vec = vdupq_n_u64(tVec[index]);

    for (long i = 0; i < N; i += 2) {
        uint64x2_t a_vec = vld1q_u64(a + i);
        uint64x2_t b_vec = vld1q_u64(b + i);

        uint64x2_t sum = vaddq_u64(a_vec, b_vec);
        uint64x2_t cmp = vcgtq_u64(sum, t_vec);
        uint64x2_t res_vec = vsubq_u64(sum, vandq_u64(cmp, t_vec));
        vst1q_u64(a + i, res_vec);
    }
}

void Context::TAddAndEqual(uint64_t *a, uint64_t *b, long t)
{
#pragma omp parallel for
    for (long i = 0; i < t; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        TiAddAndEqual(ai, bi, i);
    }
}

void Context::AddAndEqual(uint64_t *a, uint64_t *b, long l, long k)
{
#pragma omp parallel for
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);

        uint64x2_t q_vec = vdupq_n_u64(qVec[i]);
        #pragma omp parallel for
        for (long n = 0; n < N; n += 2) {
            uint64x2_t a_vec = vld1q_u64(ai + n);
            uint64x2_t b_vec = vld1q_u64(bi + n);
            uint64x2_t sum = vaddq_u64(a_vec, b_vec);
            uint64x2_t cmp = vcgtq_u64(sum, q_vec);
            uint64x2_t res_vec = vsubq_u64(sum, vandq_u64(cmp, q_vec));
            vst1q_u64(ai + n, res_vec);
        }
    }

    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        PiAddAndEqual(ai, bi, i - l);
    }
}

void Context::QiSub(uint64_t *res, uint64_t *a, uint64_t *b, long index)
{
    for (long i = 0; i < N; i++) {
        SubMod(res[i], a[i], b[i], qVec[index]);
    }
}

void Context::PiSub(uint64_t *res, uint64_t *a, uint64_t *b, long index)
{
    for (long i = 0; i < N; i++) {
        SubMod(res[i], a[i], b[i], pVec[index]);
    }
}

void Context::Sub(uint64_t *res, uint64_t *a, uint64_t *b, long l, long k)
{
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        uint64_t *resi = res + (i << logN);
        QiSub(resi, ai, bi, i);
    }

    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        uint64_t *resi = res + (i << logN);
        PiSub(resi, ai, bi, i - l);
    }
}

void Context::QiSubAndEqual(uint64_t *a, uint64_t *b, long index)
{
    uint64_t r;
    uint64_t q = qVec[index];
    for (long i = 0; i < N; i++) {
        r = b[i] % q;
        a[i] = (a[i] + q - r) % q;
    }
}

void Context::PiSubAndEqual(uint64_t *a, uint64_t *b, long index)
{
    uint64_t r;
    uint64_t pModule = pVec[index];
    for (long i = 0; i < N; i++) {
        r = b[i] % pModule;
        a[i] = (a[i] + pModule - r) % pModule;
    }
}

void Context::SubAndEqual(uint64_t *a, uint64_t *b, long l, long k)
{
#pragma omp parallel for
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        QiSubAndEqual(ai, bi, i);
    }

    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        PiSubAndEqual(ai, bi, i - l);
    }
}

void Context::TSubAndEqual(uint64_t *a, uint64_t *b, long l, long k)
{
    for (long i = 0; i < k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        PiSubAndEqual(ai, bi, i);
    }
#pragma omp parallel for
    for (long i = k; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        QiSubAndEqual(ai, bi, i - k);
    }
}

void Context::QiSub2AndEqual(uint64_t *a, uint64_t *b, long index)
{
    for (long i = 0; i < N; i++) {
        SubMod(b[i], a[i], b[i], qVec[index]);
    }
}

void Context::PiSub2AndEqual(uint64_t *a, uint64_t *b, long index)
{
    for (long i = 0; i < N; i++) {
        SubMod(b[i], a[i], b[i], pVec[index]);
    }
}

void Context::Sub2AndEqual(uint64_t *a, uint64_t *b, long l, long k)
{
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        QiSub2AndEqual(ai, bi, i);
    }

    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        PiSub2AndEqual(ai, bi, i - l);
    }
}


void Context::QiMulConstAndEqual(uint64_t *res, uint64_t cnst, long index)
{
    for (long i = 0; i < N; i++) {
        MulModBarrett(res[i], res[i], cnst, qVec[index], twoqVec[index], qrVecHigh[index], qrVecLow[index]);
    }
}

void Context::PiMulConstAndEqual(uint64_t *res, uint64_t cnst, long index)
{
    for (long i = 0; i < N; i++) {
        MulModBarrett(res[i], res[i], cnst, pVec[index], twopVec[index], prVecHigh[index], prVecLow[index]);
    }
}

void Context::MulConstAndEqual(uint64_t *res, uint64_t cnst, long l, long k)
{
    for (long i = 0; i < l; i++) {
        uint64_t *ai = res + (i << logN);
        QiMulConstAndEqual(ai, cnst, i);
    }

    for (long i = l; i < l + k; i++) {
        uint64_t *ai = res + (i << logN);
        PiMulConstAndEqual(ai, cnst, i - l);
    }
}

void Context::QiMul(uint64_t *res, uint64_t *a, uint64_t *b, long index)
{
    #pragma omp parallel for
    for (long i = 0; i < N; i++) {
        MulModBarrett(res[i], a[i], b[i], qVec[index], twoqVec[index], qrVecHigh[index], qrVecLow[index]);
    }
}

void Context::PiMul(uint64_t *res, uint64_t *a, uint64_t *b, long index)
{
    #pragma omp parallel for
    for (long i = 0; i < N; i++) {
        MulModBarrett(res[i], a[i], b[i], pVec[index], twopVec[index], prVecHigh[index], prVecLow[index]);
    }
}


void Context::Mul(uint64_t *res, uint64_t *a, uint64_t *b, long l, long k)
{
#pragma omp parallel for
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        uint64_t *resi = res + (i << logN);
        QiMul(resi, ai, bi, i);
    }
    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        uint64_t *resi = res + (i << logN);
        PiMul(resi, ai, bi, i - l);
    }
}

void Context::QiMulAndEqual(uint64_t *a, uint64_t *b, long index)
{
    uint64x2_t q_vec = vdupq_n_u64(qVec[index]);
    uint64x2_t twoq_vec = vdupq_n_u64(twoqVec[index]);
    uint64x2_t qr_high_vec = vdupq_n_u64(qrVecHigh[index]);
    uint64x2_t qr_low_vec = vdupq_n_u64(qrVecLow[index]);

    for (long i = 0; i < N; i += 2) {
        uint64x2_t a_vec = vld1q_u64(a + i);
        uint64x2_t b_vec = vld1q_u64(b + i);
        uint64x2_t res_vec;

        MulModBarrettNeon(&res_vec, a_vec, b_vec, q_vec, twoq_vec, qr_high_vec, qr_low_vec);
        vst1q_u64(a + i, res_vec);
    }
}

void Context::PiMulAndEqual(uint64_t *a, uint64_t *b, long index)
{
    uint64x2_t p_vec = vdupq_n_u64(pVec[index]);
    uint64x2_t twop_vec = vdupq_n_u64(twopVec[index]);
    uint64x2_t pr_high_vec = vdupq_n_u64(prVecHigh[index]);
    uint64x2_t pr_low_vec = vdupq_n_u64(prVecLow[index]);

    for (long i = 0; i < N; i += 2) {
        uint64x2_t a_vec = vld1q_u64(a + i);
        uint64x2_t b_vec = vld1q_u64(b + i);
        uint64x2_t res_vec;

        MulModBarrettNeon(&res_vec, a_vec, b_vec, p_vec, twop_vec, pr_high_vec, pr_low_vec);
        vst1q_u64(a + i, res_vec);
    }
}

void Context::MulAndEqual(uint64_t *a, uint64_t *b, long l, long k)
{
#pragma omp parallel for
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        QiMulAndEqual(ai, bi, i);
    }

    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *bi = b + (i << logN);
        PiMulAndEqual(ai, bi, i - l);
    }
}

void Context::QiSquare(uint64_t *res, uint64_t *a, long index)
{
    for (long i = 0; i < N; i++) {
        MulModBarrett(res[i], a[i], a[i], qVec[index], twoqVec[index], qrVecHigh[index], qrVecLow[index]);
    }
}

void Context::PiSquare(uint64_t *res, uint64_t *a, long index)
{
    for (long i = 0; i < N; i++) {
        MulModBarrett(res[i], a[i], a[i], pVec[index], twopVec[index], prVecHigh[index], prVecLow[index]);
    }
}

void Context::Square(uint64_t *res, uint64_t *a, long l, long k)
{
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *resi = res + (i << logN);
        QiSquare(resi, ai, i);
    }

    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        uint64_t *resi = res + (i << logN);
        PiSquare(resi, ai, i - l);
    }
}


void Context::SquareAndEqual(uint64_t *a, long l, long k)
{
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        QiSquare(ai, ai, i);
    }

    for (long i = l; i < l + k; i++) {
        uint64_t *ai = a + (i << logN);
        PiSquare(ai, ai, i - l);
    }
}

void Context::EvalAndEqual(uint64_t *a, long l)
{
    INTTAndEqual(a, l);
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        for (long n = 0; n < N; n++) {
            MulModBarrett(ai[n], ai[n], PModq[i], qVec[i], twoqVec[i], qrVecHigh[i], qrVecLow[i]);
        }
    }
    NTTAndEqual(a, l);
}

void Context::Raise(uint64_t *res, uint64_t *a, long l) {}

void Context::RaiseAndEqual(uint64_t *&a, long l)
{
    std::copy(a, a + (l << logN), raRaiseAndEqual);

    INTTAndEqual(a, l);

    for (long i = 0; i < l; i++) {
        uint64_t *tmp3i = tmp3 + (i << logN);
        uint64_t *ai = a + (i << logN);
        for (long n = 0; n < N; n++) {
            MulModBarrett(tmp3i[n], ai[n], qHatInvModq[l - 1][i], qVec[i], twoqVec[i], qrVecHigh[i], qrVecLow[i]);
        }
    }
    for (long k = 0; k < K; k++) {
        uint64_t *rak = raRaiseAndEqual + ((l + k) << logN);
        for (long n = 0; n < N; n++) {
            uint64_t tt = tmp3[n];
            uint128_t sum = static_cast<uint128_t>(tt) * qHatModp[l - 1][0][k];
            for (long i = 1; i < l; i++) {
                tt = tmp3[n + (i << logN)];
                sum += static_cast<uint128_t>(tt) * qHatModp[l - 1][i][k];
            }
            ModBarrett(rak[n], sum, pVec[k], prVecHigh[k], prVecLow[k]);
        }
    }
    NTTAndEqual(raRaiseAndEqual + (l << logN), 0, K);

    delete[] a;
    a = raRaiseAndEqual;
}

void Context::Back(uint64_t *res, uint64_t *a, long l)
{
    std::copy(a, a + ((l + K) << logN), tmpCode);
    INTTAndEqual(tmpCode, l, K);
    for (long k = 0; k < K; k++) {
        uint64_t *tmpk = tmpCode + ((l + k) << logN);
        uint64_t *tmp3k = tmp3 + (k << logN);
        for (long n = 0; n < N; n++) {
            MulModBarrett(tmp3k[n], tmpk[n], pHatInvModp[k], pVec[k], twopVec[k], prVecHigh[k], prVecLow[k]);
        }
    }
    for (long i = 0; i < l; i++) {
        uint64_t *resi = res + (i << logN);
        uint64_t *tmpi = tmpCode + (i << logN);

        for (long n = 0; n < N; n++) {
            uint64_t tt = tmp3[n];
            uint128_t sum = static_cast<uint128_t>(tt) * pHatModq[0][i];
            for (long k = 1; k < K; k++) {
                tt = tmp3[n + (k << logN)];
                sum += static_cast<uint128_t>(tt) * pHatModq[k][i];
            }
            ModBarrett(resi[n], sum, qVec[i], qrVecHigh[i], qrVecLow[i]);
            SubMod(resi[n], tmpi[n], resi[n], qVec[i]);
            MulModBarrett(resi[n], resi[n], PInvModq[i], qVec[i], twoqVec[i], qrVecHigh[i], qrVecLow[i]);
        }
    }

    NTTAndEqual(res, l);
}

void Context::BackAndEqual(uint64_t *&a, long l)
{
    INTTAndEqual(a, l, K);


    for (long k = 0; k < K; k++) {
        uint64_t *tmp3k = tmp3 + (k << logN);
        uint64_t *ak = a + ((k + l) << logN);
        for (long n = 0; n < N; n++) {
            MulModBarrett(tmp3k[n], ak[n], pHatInvModp[k], pVec[k], twopVec[k], prVecHigh[k], prVecLow[k]);
        }
    }

    for (long i = 0; i < l; i++) {
        uint64_t *rai = raBackAndEqual + (i << logN);
        uint64_t *ai = a + (i << logN);
        for (long n = 0; n < N; n++) {
            uint64_t tt = tmp3[n];
            uint128_t sum = static_cast<uint128_t>(tt) * pHatModq[0][i];
            for (long k = 1; k < K; k++) {
                tt = tmp3[n + (k << logN)];
                sum += static_cast<uint128_t>(tt) * pHatModq[k][i];
            }
            ModBarrett(rai[n], sum, qVec[i], qrVecHigh[i], qrVecLow[i]);
            SubMod(rai[n], ai[n], rai[n], qVec[i]);
            MulModBarrett(rai[n], rai[n], PInvModq[i], qVec[i], twoqVec[i], qrVecHigh[i], qrVecLow[i]);
        }
    }
    NTTAndEqual(raBackAndEqual, l);
    delete[] a;
    a = raBackAndEqual;
}

void Context::ReScale(uint64_t *res, uint64_t *a, long l)
{
    uint64_t *al = a + ((l - 1) << logN);
    for (long i = 0; i < l - 1; i++) {
        uint64_t *rai = raReScale + (i << logN);
        uint64_t *ai = a + (i << logN);

        std::copy(al, al + N, rai);
        for (long n = 0; n < N; n++) {
            ModBarrett(rai[n], al[n], qVec[i], qrVecHigh[i], qrVecLow[i]);
        }
        QiNTTAndEqual(rai, i);

        for (long n = 0; n < N; n++) {
            SubMod(rai[n], ai[n], rai[n], qVec[i]);
            MulModBarrett(rai[n], rai[n], qInvModq[l - 1][i], qVec[i], twoqVec[i], qrVecHigh[i], qrVecLow[i]);
        }
    }

    std::copy(raReScale, raReScale + ((l - 1) << logN), res);
}

void Context::ReScaleAndEqual(uint64_t *a, long l)
{
    uint64_t *ra = raReScaleAndEqual;
    uint64_t *al = a + ((l - 1) << logN);
    QiINTTAndEqual(al, l - 1);
    #pragma omp parallel for
    for (long i = 0; i < l - 1; i++) {
        uint64_t *rai = ra + (i << logN);
        uint64_t *ai = a + (i << logN);
        std::copy(al, al + N, rai);
        #pragma omp parallel for
        for (long n = 0; n < N; n++) {
            ModBarrett(rai[n], al[n], qVec[i], qrVecHigh[i], qrVecLow[i]);
        }
        QiNTTAndEqual(rai, i);

        #pragma omp parallel for
        for (long n = 0; n < N; n++) {
            SubMod(rai[n], ai[n], rai[n], qVec[i]);
            MulModBarrett(rai[n], rai[n], qInvModq[l - 1][i], qVec[i], twoqVec[i], qrVecHigh[i], qrVecLow[i]);
        }
    }
    std::copy(ra, ra + ((l - 1) << logN), a);
}

void Context::ModDownAndEqual(uint64_t *&a, long l, long dl)
{
    uint64_t *ra = new uint64_t[(l - dl) << logN];
    std::copy(a, a + ((l - dl) << logN), ra);
    delete[] a;
    a = ra;
}

uint64_t *Context::ModDown(uint64_t *a, long l, long dl)
{
    uint64_t *ra = new uint64_t[(l - dl) << logN];
    std::copy(a, a + ((l - dl) << logN), ra);
    return ra;
}

void Context::LeftRot(uint64_t *res, uint64_t *a, long l, long rotSlots)
{
    std::copy(a, a + (l << logN), tmpLeftRot);
    INTTAndEqual(tmpLeftRot, l);
    long pow = rotGroup[rotSlots];
    for (long i = 0; i < l; i++) {
        uint64_t *resi = res + (i << logN);
        uint64_t *tmpi = tmpLeftRot + (i << logN);
        for (long n = 0; n < N; n++) {
            long nPow = n * pow;
            long shift = nPow % M;
            if (shift < N) {
                resi[shift] = tmpi[n];
            } else {
                resi[shift - N] = qVec[i] - tmpi[n];
            }
        }
    }
    NTTAndEqual(res, l);
}

void Context::LeftRotAndEqual(uint64_t *a, long l, long rotSlots)
{
    std::copy(a, a + (l << logN), tmpCode);
    long idx = rotSlots % Nh;
    for (long n = 0; n < N; n++) {
        uint32_t reversed = BitReverse(static_cast<uint32_t>(n)) >> (32 - logN);
        uint64_t index_raw = rotGroup[idx] * (2 * reversed + 1);
        index_raw &= M - 1;
        uint32_t index = BitReverse((static_cast<uint32_t>(index_raw) - 1) >> 1) >> (32 - logN);
        for (long i = 0; i < l; i++) {
            a[n + (i << logN)] = tmpCode[index + (i << logN)];
        }
    }
}

void Context::Conjugate(uint64_t *res, uint64_t *a, long l)
{
    for (long i = 0; i < l; i++) {
        uint64_t *resi = res + (i << logN);
        uint64_t *ai = a + (i << logN);
        for (int n = 0; n < N; n++) {
            resi[n] = ai[N - 1 - n];
        }
    }
}

void Context::ConjugateAndEqual(uint64_t *a, long l)
{
    for (long i = 0; i < l; i++) {
        uint64_t *ai = a + (i << logN);
        for (int n = 0; n < N; n++) {
            swap(ai[n], ai[N - 1 - n]);
        }
    }
}

/*
 * @brief ：生成服从高斯分布的多项式系数样本
 *
 * @param res： 输出数组，用于存储生成的高斯分布样本。
 * @param l： 处理的模数层级数量，前 `l` 层使用标准模数向量 `qVec`。
 * @param k： 额外层级数量，后 `k` 层使用额外模数向量 `pVec`。
 */
void Context::SampleGauss(uint64_t *res, long l, long k)
{
    // 使用 C++11 <random> 库的随机数生成器和分布
    // 创建静态随机数生成器，确保只初始化一次
    static std::mt19937_64 gen(std::random_device{}());

    // 定义均匀分布，范围 (0, 1)
    static std::uniform_real_distribution<double> dis(0.0, 1.0);

    for (long i = 0; i < N; i += 2) {
        // 生成两个独立的均匀分布随机数，范围 (0, 1)
        double r1 = dis(gen);
        double r2 = dis(gen);

        // 使用 Box-Muller 变换将均匀分布转换为标准正态分布
        double theta = 2.0 * M_PI * r1;
        double rr = sqrt(-2.0 * log(r2)) * sigma;

        // 生成两个高斯随机数，并四舍五入到最近的整数
        long g1 = static_cast<long>(std::floor(rr * std::cos(theta) + 0.5));
        long g2 = static_cast<long>(std::floor(rr * std::sin(theta) + 0.5));

        // 将生成的高斯随机数应用到结果数组
        for (long j = 0; j < l; j++) {
            uint64_t *resj = res + (j << logN);
            resj[i] = g1 >= 0 ? static_cast<uint64_t>(g1) : qVec[j] + g1;
            resj[i + 1] = g2 >= 0 ? static_cast<uint64_t>(g2) : qVec[j] + g2;
        }
        for (long j = 0; j < k; j++) {
            uint64_t *resj = res + ((j + l) << logN);
            resj[i] = g1 >= 0 ? static_cast<uint64_t>(g1) : pVec[j] + g1;
            resj[i + 1] = g2 >= 0 ? static_cast<uint64_t>(g2) : pVec[j] + g2;
        }
    }
}

/*
 * @brief ：生成服从 {0, 1, -1} 分布的多项式系数样本
 *
 * @param res： 输出数组，用于存储生成的 {0, +1, -1} 分布样本。
 * @param s： 未使用的参数（预留接口一致性，可能为以后扩展所用）。
 * @param l： 处理的模数层级数量，前 `l` 层使用标准模数向量 `qVec`。
 * @param k： 额外层级数量，后 `k` 层使用额外模数向量 `pVec`。
 */
void Context::SampleZO(uint64_t *res, long l, long k)
{
    // 使用 C++11 <random> 库的随机数生成器和分布
    // 创建静态随机数生成器，确保只初始化一次
    static std::mt19937_64 gen(std::random_device{}());

    // 定义分布
    // 对于 zo = {0, 1, -1}，每个值的概率均为1/3
    std::uniform_int_distribution<int> dist_zo(0, 2); // 0,1,2 分别映射到 0,1,-1

    for (long i = 0; i < N; i++) {
        int rand_val = dist_zo(gen);
        long zo;
        switch (rand_val) {
            case 0:
                zo = 0;
                break;
            case 1:
                zo = 1;
                break;
            case 2:
                zo = -1;
                break;
            default:
                zo = 0; // 理论上不会到达
        }

        if (std::abs(zo) > 1)
            std::cerr << "error in TUG !!!!!!!\n";

        for (long j = 0; j < l; j++) {
            uint64_t *resj = res + (j << logN);
            resj[i] = zo >= 0 ? static_cast<uint64_t>(zo) : qVec[j] + zo;
        }
        for (long j = 0; j < k; j++) {
            uint64_t *resj = res + ((j + l) << logN);
            resj[i] = zo >= 0 ? static_cast<uint64_t>(zo) : pVec[j] + zo;
        }
    }
}

/*
 * @brief ：生成均匀分布的随机多项式系数样本。
 *
 * @param res： 输出数组，用于存储生成的均匀分布样本。
 * @param l： 处理的模数层级数量，前 `l` 层使用标准模数向量 `qVec`。
 * @param k： 额外层级数量，后 `k` 层使用额外模数向量 `pVec`。
 */
void Context::SampleUniform(uint64_t *res, long l, long k)
{
    // 使用 C++11 <random> 库的随机数生成器和分布
    // 创建静态随机数生成器，确保只初始化一次
    static std::mt19937_64 gen(std::random_device{}());

    // 定义均匀实数分布，范围 [0.0, 1.0)
    std::uniform_real_distribution<double> dis(0.0, 1.0);

    // 生成均匀分布的随机数并赋值
    for (long j = 0; j < l; j++) {
        uint64_t *resj = res + (j << logN);
        for (long n = 0; n < N; n++) {
            double rand_real = dis(gen);
            resj[n] = static_cast<uint64_t>(std::floor(rand_real * qVec[j]));
        }
    }
    for (long j = 0; j < k; j++) {
        uint64_t *resj = res + ((j + l) << logN);
        for (long n = 0; n < N; n++) {
            double rand_real = dis(gen);
            resj[n] = static_cast<uint64_t>(std::floor(rand_real * pVec[j]));
        }
    }
}

void Context::SampleHWT(uint64_t *res, long l, long k)
{
    // 使用 C++11 <random> 库的随机数生成器和分布
    // 创建静态随机数生成器，确保只初始化一次
    static std::mt19937_64 gen(std::random_device{}());

    // 定义均匀实数分布，范围 [0.0, 1.0)
    std::uniform_real_distribution<double> dis_real(0.0, 1.0);

    // 定义均匀整数分布，范围 [0, N-1]
    std::uniform_int_distribution<long> dis_int(0, N - 1);

    // 定义均匀整数分布，范围 [0, 1] 映射到 {1, -1}
    std::uniform_int_distribution<int> dis_sign(0, 1);

    long idx = 0;
    while (idx < h) {
        long i = dis_int(gen);

        if (res[i] == 0) {
            int sign = dis_sign(gen);
            long hwt = sign ? 1 : -1;

            for (long j = 0; j < l; j++) {
                uint64_t *resj = res + (j << logN);
                resj[i] = hwt >= 0 ? static_cast<uint64_t>(hwt) : qVec[j] + hwt;
            }
            for (long j = 0; j < k; j++) {
                uint64_t *resj = res + ((j + l) << logN);
                resj[i] = hwt >= 0 ? static_cast<uint64_t>(hwt) : pVec[j] + hwt;
            }
            idx++;
        }
    }
}

/*
 * @brief ：计算 (a^k) % mod 的快速幂结果，使用二进制快速幂算法。
 *
 * @param a： 底数，参与指数幂运算的数值。
 * @param k： 指数，表示幂次。
 * @param mod： 模数，用于取模运算，确保结果在 [0, mod) 范围内。
 *
 * @return uint64_t: (a^k) % mod 的计算结果。
 */
uint64_t Context::FastPow(uint64_t a, uint64_t k, uint64_t mod)
{
    uint64_t base = 1;
    __uint128_t b = a;
    while (k) {
        if (k & 1) {
            if (mod == 0) {
                throw std::runtime_error("Divisor mod by zero is not allowed.");
            }
            base = (base * b) % mod;
        }
        b = (b * b) % mod;
        k >>= 1;
    }
    return base % mod;
}


uint64_t Context::BitReverse64(uint64_t index, uint64_t bitLen)
{
    uint64_t reversed = 0;
    for (uint64_t i = 0; i < bitLen; i++) {
        reversed = (reversed << 1) | (index & 1);
        index >>= 1;
    }
    return reversed;
}

/*
 * @brief ：生成 NTT 索引重排数组，基于输入的伽罗瓦元素进行重排。
 *
 * @param galEl ：输入的伽罗瓦元素，用于计算索引重排的映射关系。
 * @param ntt_index ：输出数组，用于存储重排后的 NTT 索引。
 */
void Context::PermuteNTTIndex(uint64_t galEl, uint64_t *ntt_index)
{
    uint64_t logNthRoot = log2(M) - 1;
    uint64_t mask = M - 1;
    #pragma omp parallel for
    for (int i = 0; i < N; i++) {
        uint64_t tmp1 = 2 * BitReverse64(i, logNthRoot) + 1;
        uint64_t tmp2 = ((galEl * tmp1 & mask) - 1) >> 1;
        ntt_index[i] = BitReverse64(tmp2, logNthRoot);
    }
}

/*
 * @brief ：使用给定的重排索引数组，对输入多项式数组执行重排操作。
 *
 * @param polOut： 输出数组，用于存储经过重排后的多项式系数。
 * @param polIn： 输入数组，包含待重排的多项式系数。
 * @param index： 重排索引数组，定义多项式系数的重排顺序。
 * @param level： 层级数量，表示要处理的多项式分块数量。
 */
void Context::PermuteNTTWithIndex(uint64_t *polOut, uint64_t *polIn, uint64_t *index, int level)
{
    for (long l = 0; l < level; l++) {
        uint64_t *polIni = polIn + (l << logN);
        uint64_t *polOuti = polOut + (l << logN);
        IPermuteNTTWithIndex(polOuti, polIni, index);
    }
}

/*
 * @brief ：使用指定的重排索引数组，对输入多项式数组进行重新排列。
 *
 * @param polOut： 输出数组，用于存储重排后的多项式系数。
 * @param polIn： 输入数组，包含待重排的多项式系数。
 * @param index： 重排索引数组，指定多项式系数的重排顺序。
 */
void Context::IPermuteNTTWithIndex(uint64_t *polOut, uint64_t *polIn, uint64_t *index)
{
    #pragma omp parallel for
    for (int i = 0; i < N; i++) {
        polOut[i] = polIn[index[i]];
    }
}


/*
 * @brief ：执行基数转换操作，将输入多项式数组从模数 PQ 转换到目标模数 T 上
 *
 * @param a ：输入数组，包含待转换的多项式系数。
 * @return vector<vector<uint64_t>> ：返回一个二维向量，存储转换后的多项式系数（每个向量对应一个目标模数）。
 */
vector<vector<uint64_t>> Context::BaseConvertPQ2T(uint64_t *a)
{
    vector<vector<uint64_t>> res;
    for (int i = 0; i < blockLen; i++) {
        vector<uint64_t> convT;
        convT.resize(T << logN);
        res.emplace_back(convT);
    }
    #pragma omp parallel for collapse(2)
    for (int i = 0; i < blockLen; i++) {
        for (int n = 0; n < N; n++) {
            int index = i * gamma;
            int length = gamma;
            if (i == blockLen - 1) {
                length = K + L - index;
            }
            float vv = 0;
            uint64_t ra[BASECONVERT];
            for (int j = 0; j < length; j++) {
                uint64_t *aij = a + ((index + j) << logN);
                    MulModShoup(ra[j], aij[n], PQSlicesHatInvModPQ[i][j], PQSlicesHatInvModPQShoup[i][j],
                                PQ[index + j]);
                    vv += ra[j] / float(PQ[index + j]);
            }
            for (int j = 0; j < T; j++) {
                int offset = j << logN;
                uint128_t sum = 0;
                for (int m = 0; m < length; m++) {
                    uint64_t tt = ra[m];
                    sum += static_cast<uint128_t>(tt) * PQSlicesHatModt[i][m][j];
                }
                uint64_t PQMul = int(vv) * PQSlicesModt[i][j];
                ModBarrett(res[i][offset + n], sum - PQMul, tVec[j], trVecHigh[j], trVecLow[j]);
            }
        }
    }

    for (int i = 0; i < blockLen; i++) {
        TNTTAndEqual(res[i].data(), T);
    }
   return res;
}


vector<uint64_t *> Context::BaseConvertQl2T(uint64_t *a, long l)
{
    float numSlices = std::ceil(l / float(K));
    vector<uint64_t *> res;

    for (int i = 0; i < numSlices; i++) {
        int offset = i * K;
        int length = K;
        if (i == numSlices - 1) {
            length = l - offset;
        }
        uint64_t *ai = a + (offset << logN);
        res.emplace_back(SlicesBaseConvertQl2T(ai, length, i, offset));
    }
    return res;
}

// /*
//  * @brief ：将输入多项式数组从模数 Q 分组转换到目标模数 T
//  *
//  * @param a： 输入数组，包含模数集合 Q 下的多项式系数。
//  * @param l： 模数 Q 的层级数量。
//  * @param res： 输出结果向量，存储每个分块转换后的结果数组指针。
//  */

// /*
//  * @brief 执行单个分块的模数集合 Ql 到目标模数集合 T 的基数转换。
//  *
//  * @param a： 输入数组，包含模数集合 Ql 下的多项式系数。
//  * @param length： 当前分块的长度，表示模数层级数量。
//  * @param slices： 当前分块的索引。
//  * @param offset： 当前分块在整体模数集合中的起始偏移。
//  *
//  * @return uint64_t*： 返回转换后的结果数组指针，存储模数集合 T 下的多项式系数。
//  */
void Context::BaseConvertQl2T(uint64_t *a, long l, vector<uint64_t *> &res)
{
    res.clear();
    float numSlices = std::ceil(l / float(K));
    for (int i = 0; i < numSlices; i++) {
        int offset = i * K;
        int length = K;
        if (i == numSlices - 1) {
            length = l - offset;
        }
        uint64_t *ai = a + (offset << logN);
        res.emplace_back(SlicesBaseConvertQl2T(ai, length, i, offset));
    }
}

// /*
//  * @brief 执行单个分块的模数集合 Ql 到目标模数集合 T 的基数转换。
//  *
//  * @param a： 输入数组，包含模数集合 Ql 下的多项式系数。
//  * @param length： 当前分块的长度，表示模数层级数量。
//  * @param slices： 当前分块的索引。
//  * @param offset： 当前分块在整体模数集合中的起始偏移。
//  *
//  * @return uint64_t*： 返回转换后的结果数组指针，存储模数集合 T 下的多项式系数。
//  */
uint64_t *Context::SlicesBaseConvertQl2T(uint64_t *a, long length, long slices, long offset)
{
    uint64_t *convT = new uint64_t[T << logN];
#pragma omp parallel for
    for (int n = 0; n < N; n++) {
        uint64_t ra[BASECONVERT];
        float vv = 0;
        for (int j = 0; j < length; j++) {
            uint64_t *aj = a + (j << logN);
            MulModBarrett(ra[j], aj[n], qSlicesHatInvModq[slices][length - 1][j], qVec[offset + j], twoqVec[offset + j], qrVecHigh[offset + j], qrVecLow[offset + j]);
            vv += ra[j] * qVecInv[offset + j];;
        }
        for (int j = 0; j < T; j++) {
            uint64_t jN = j << logN;
            uint128_t sum = 0;
            for (int m = 0; m < length; m++) {
                __uint128_t tt = ra[m];
                sum += tt * qSlicesHatModt[slices][length - 1][m][j];
            }
            uint64_t QMul = static_cast<int>(vv) * QSlicesModt[slices][length - 1][j];
            ModBarrett(convT[jN + n], sum - QMul, tVec[j], trVecHigh[j], trVecLow[j]);
        }
    }
    return convT;
}


// /*
//  * @brief ：将输入多项式系数从模数 T 转换到目标模数 PQ。
//  *
//  * @param res： 输出数组，存储转换后的多项式系数。
//  * @param a： 输入数组指针，包含模数 T 下的多项式数据。
//  * @param l： 模数 PQl 的层级数量，用于计算转换的分块长度。
//  */
void Context::BaseConvertT2PQl(uint64_t *res, uint64_t **a, long l)
{
    uint64_t TMul;
    int numSlices = std::ceil((K + l) / float(gamma));
#pragma omp parallel for collapse(2)
    for (int d = 0; d < numSlices; d++) {
        for (int n = 0; n < N; n++) {
            int offset = d * gamma;
            int length = gamma;
            if (d == numSlices - 1) {
                length = K + l - offset;
            }

            uint64_t ra[BASECONVERT];
            float vv = 0;
            for (int i = 0; i < T; i++) {
                uint64_t *ai = a[d] + (i << logN);
                MulModShoup(ra[i], ai[n], tHatInvModt[i], tHatInvModtShoup[i], tVec[i]);
                vv += ra[i] / float(tVec[i]);
            }

            for (int j = 0; j < length; j++) {
                uint64_t djN = (offset + j) << logN;
                uint128_t sum = 0;
                for (int m = 0; m < T; m++) {
                    uint64_t tt = ra[m];
                    sum += static_cast<uint128_t>(tt) * tHatModPQ[m][offset + j];
                }
                TMul = int(vv) * TModPQ[offset + j];
                ModBarrett(res[djN + n], sum - TMul, PQ[offset + j], pqrVecHigh[offset + j],
                    pqrVecLow[offset + j]);
            }
        }
    }
}

// /*
//  * @brief ：将输入多项式系数从模数 T 转换到目标模数 PQl。
//  *
//  * @param a： 输入数组指针，包含模数集合 T 下的多项式数据。
//  * @param l： 模数集合 PQl 的层级数量，用于计算转换的分块长度。
//  * @return uint64_t* ：输出数组，存储转换后的多项式系数。
//  */
void Context::BaseConvertPQl2Ql(uint64_t *res, uint64_t *a, long l)
{
    if( K != 1) {
        #pragma omp parallel for collapse(2)
        for (int i = 0; i < K; i++) {
            for (int n = 0; n < N; n += 2) {
                uint64_t *ai = a + (i << logN);
                uint64x2_t ai_vec = vld1q_u64(ai + n);
                uint64x2_t pHatInvModp_vec = vdupq_n_u64(pHatInvModp[i]);
                uint64x2_t pHatInvModpShoup_vec = vdupq_n_u64(pHatInvModpShoup[i]);
                uint64x2_t pVec_vec = vdupq_n_u64(pVec[i]);

                ShoupMulNEON(&ai_vec, ai_vec, pHatInvModp_vec, pHatInvModpShoup_vec, pVec_vec);
                vst1q_u64(ai + n, ai_vec);
            }
        }
    }
    
#pragma omp parallel for collapse(2)
    for (int i = 0; i < l; i++) {
        for (int n = 0; n < N; n += 2) {
            int offset = i << logN;
            uint64x2_t baseConvPQ2Q_vec = vld1q_u64(a + n);
            __uint128_t sum_low = static_cast<__uint128_t>(vgetq_lane_u64(baseConvPQ2Q_vec, 0)) * pHatModq[0][i];
            __uint128_t sum_high = static_cast<__uint128_t>(vgetq_lane_u64(baseConvPQ2Q_vec, 1)) * pHatModq[0][i];

            for (int m = 1; m < K; m++) {
                baseConvPQ2Q_vec = vld1q_u64(a + n + (m << logN));
                sum_low += static_cast<__uint128_t>(vgetq_lane_u64(baseConvPQ2Q_vec, 0)) * pHatModq[m][i];
                sum_high += static_cast<__uint128_t>(vgetq_lane_u64(baseConvPQ2Q_vec, 1)) * pHatModq[m][i];
            }

            uint64_t result_low, result_high;
            ModBarrett(result_low, sum_low, qVec[i], qrVecHigh[i], qrVecLow[i]);
            ModBarrett(result_high, sum_high, qVec[i], qrVecHigh[i], qrVecLow[i]);

            res[offset + n] = result_low;
            res[offset + n + 1] = result_high;
        }
    }

    uint64_t *aK = a + (K << logN);

#pragma omp parallel for collapse(2)
    for (int i = 0; i < l; i++) {
        for (int n = 0; n < N; n += 2) {
            int offset = i << logN;
            uint64x2_t aK_vec = vld1q_u64(aK + offset + n);
            uint64x2_t res_vec = vld1q_u64(res + offset + n);

            int64x2_t mid_vec = vreinterpretq_s64_u64(vsubq_u64(aK_vec, res_vec));
            while (vgetq_lane_s64(mid_vec, 0) < 0) {
                mid_vec = vaddq_s64(mid_vec, vdupq_n_s64(qVec[i]));
            }
            while (vgetq_lane_s64(mid_vec, 1) < 0) {
                mid_vec = vaddq_s64(mid_vec, vdupq_n_s64(qVec[i]));
            }

            uint64x2_t mid_u64_vec = vreinterpretq_u64_s64(mid_vec);
            uint64x2_t tmp_res_vec;
            uint64x2_t PInvModq_vec = vdupq_n_u64(PInvModq[i]);
            uint64x2_t PInvModqShoup_vec = vdupq_n_u64(PInvModqShoup[i]);
            uint64x2_t qVec_vec = vdupq_n_u64(qVec[i]);

            ShoupMulNEON(&tmp_res_vec, mid_u64_vec, PInvModq_vec, PInvModqShoup_vec, qVec_vec);
            vst1q_u64(res + offset + n, tmp_res_vec);
        }
    }
}
